disp ( 1 Program Plasma6 .m ? ) ; 

%; plasma6.m - program to compute current required for fusion in bottle 



compute the ring equation 

- int(dB) = int(i ds x r / r A 2) 



% 
% 

% 
% 
% 
% 
% 
% 
% 

fid 



- Therefore 

- int(dBx) = 

- int(dBy) = 

- int(dBz) = 



int(dsy*rz - ry*dsz/ 
int(dsx*rz - rx*dsz/ 



s 2) 
s 2) 
^2) 



int(dsx*ry - rx*dsy/ ^ 
convert each of these into an f (x, y, z, rl, pxl, pyl, pzl, thl, th2) 
polynomial matrix: each term A 4 A 3 A 2 A l - 0.5 7.0 35.0 150.0 
B vector: i ds x r/ r A 2 for 64 different points for each 3 var. 
int (dBx, dBy, dBz) = (expressed in terms of vectors) 
fopen( 1 g: /drawl /plasma6 . txt 1 , f w f ) ; 
fprintf (f id, 1 \nPlasma6.m - OutputXn 1 ); 
syms xl yl zl xO yO zO thl rl atmp; 

Bx By Bz B uO; 



dBx dBy dBz 
1/(10 A 7); 
rl*cos (thl) + xO; 
rl*sin(thl) + yO; 
zO; 

x2 y2 z2 rol th2 thel; 
= sqrt((xl - xO) A 2 + (yl - yO) 
th2 = atan((yl - y0)/(xl - xO) ) ; 
thel = acos((zl - z0)/(th2)); 
rol*cos (th2) *sin(thel) ; 
rol*sin(th2) * sin (thel) ; 
rol*cos ( thel ) ; 
rll rlx" r±y rlz; 
pzl; 
x2); 
y2); 
z2) ; 



syms 
uO = 
xl = 
yl = 
zl = 
syms 
rol 



x2 = 
y2 = 
z2 = 
syms 
syms 



s 2 + (zl - zO) A 2); 



pxl pyl 
lx = (pxl - 
ly = (pyl - 
rlz = (pzl - 
rll = sqrt(rlx A 2 + 
syms coefrl coefr2 
rl = 1.0; 



rly A 2 + rlz A 2) ; 

coefr3 ccre-fr4 coefrS arl brl rllap; 



thl 
pxl 
pyl 
pzl 
arl 
brl 
pxl 
arl 
brl 
pxl 
arl 
brl 
pxl 
arl 
brl 



0.0; 
2.0; 
0.0; 
0.0; 

[rlx A 4 rlx A 3 
[eval(rll) ;] ; 
8; 



rlx A 2 rlx A l; ] ; 



[arl; 
[brl; 
16; 
[arl; 
[brl; 
4.8; 
[arl; 
[brl; 



rlx A 4 rlx A 3 
eval(rll);] ( 

rlx A 4 rlx A 3 
eval(rll) ;], 

rlx A 4 rlx A 3 
eval(rll) ;] 



rlx A 2 rlx A l;] 



rlx A 2 rlx A l;] 



rlx A 2 rlx A l;] 



coefrl = pinv(eval (arl) ) *brl; 



disp( 'coefrl 
pxl 
pyl 
pzl 
rl 
thl 

rl 

irl 
^hl 
arl 



brl 
1;] 



= 2.0; 
= 0.0; 
= 0.0; 
= 1.0; 

3.14159/2; 
[thl A 4 thl 
[eval (rll) 
3.14159; 

[arl; thl A 4 thl A 3 thl A 2 thl A l;]; 
[brl; eval (rll) - coefrl (1) *rl A 4 



N 3 thl A 2 thl A l;] ; 
- coefrl (1) *rl A 4 



- coefrl (2) *rl A 3 - coefrl (3) *rl A 2 - coefrl (4 ) *rl A l; ] ; 



~ coefrl (2) *rl A 3 - coefrl (3) *rl A 2 - coefrl (4 ) *rl A 



thl = 3*3.1415-9/2; 

arl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 

brl = [brl; eval(rll) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coefrl (3) *rl A 2 -r coefrl (4) *rl A 
1;]; 

•hi = 2*3.14159; 
.,rl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 

brl = [brl; evai(rll) - coefrl (1) *rl A 4 - coefrl (2 ) *rl A 3 - coefrl (3) *rl A 2 - coef rl ( 4 ) *rl A 
l;]; 

coefr2 = pinv (arl) *real (brl) ; 

disp( , coefr2. . . ■ ) ; 

rl = 1.0; 

pyl = 0.0; 

pzl = 0.0; 

thl = 0.0; 

pxl = 2; 

arl - [pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = eval(rll) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coefrl (3) *rl A 2 - coef rl (4 ) *rl A l; 
brl = [atmp - coef r2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; ] ; 
pxl = 8; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = eval(rll) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coefrl (3) *rl A 2 - coef rl (4 ) *rl A l; 
brl = [brl; atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coefr2 (3) *thl A 2 - coef r2 (4 ) *thl A 
l;]; 

pxl = 18; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = eval(rll) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4) *rl A l; 
brl = [brl; atmp - coef r2 ( 1 ) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A 
l;]; 

pxl = 48; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = eval(rll) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3 ) *rl A 2 - coef rl (4 ) *rl A l; 
brl = [brl; atmp - coef r2 ( 1 ) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 ( 4 ) *thl A 
;]; 

oefr3 = pinv(arl) *real (brl)f 
disp ( 1 coefr3. . . 1 ) ; 
rl = 1.0; 
pxl =0.0; 
pzl = 0.0; 
thl = 0.0; 
pyl - 2; 

arl = [pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = eval(rll) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
brl = [atmp - coefr3 (1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 ( 4 ) *pxl A l; ] ; 
pyl = 8; 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = eval(rll) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4) *rl A l; 
atmp - atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4) *thl A l; 
brl = [brl; atmp - coef r3 (1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A 

l;]; 

pyl = 18; 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = eval(rll) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coefrl (3) *rl A 2 - coefrl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
brl = [brl; atmp - coef r3 (1) *pxl A 4 - coefr3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 ( 4 ) *pxl A 

l;]; 

pyl = 48; 

arl = [arl; pyl A 4 pyl A 3 pyl"2 pyl A l;]; 

atmp = eval(rll) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coefrl (3) *rl A 2 - coefrl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
^Drl = [brl; atmp - coef r3 ( 1 ) *pxl A 4 - coef r3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 ( 4 ) *pxl A 

• ;]; 

coefr4 = pinv(arl) *real (brl) ; 
disp ( 1 coefr4 ...'); 
pxl = 0.0; 
pyl =0.0; 



thl = 0.0; 
rl = 1.0; 
pzl = 2; 

arl = [pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

•tmp = eval(rll) - coef rl ( 1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
tmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4 ) *thl A l; 
atmp = atmp - coefr3 (1) *pxl A 4 - coef r3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
brl = [atmp - coefr4 (1) *pyl A 4 - coef r4 (2) *pyl A 3 - coef r4 (3) *pyl A 2 - coef r4 (4 ) *pyl A l; ] ; 
pzl = 8; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = eval(rll) - coef rl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4) *thl A l; 
atmp = atmp - coefr3 (1) *pxl A 4 - coefr3 (2) *pxl A 3 - coefr3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
brl = [brl; atmp - coefr4 (1) *pyl A 4 - coef r4 (2 ) *pyl A 3 - coef r4 (3) *pyl A 2 - coef r4 (4 ) *pyl A 
1;]; 

pzl = 32; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp - eval(rll) - coef rl ( 1) *rl A 4 - coef rl (2) *rl A 3 - coefrl (3) *rl A 2 - coefrl (4) *rl A l; 
atmp « atmp - coefr2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
atmp = atmp - coefr3 (1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coefr3 (4 ) *pxl A l; 
brl « [brl; atmp - coef r4 (1) *pyl A 4 - coef r4 (2 ) *pyl A 3 - coef r4 (3) *pyl A 2 - coef r4 (4 ) *pyl A 
1;]; 

pzl = 48; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = evcri(rll) - coef rl (l)^rl A 4 - coefrl (2) *rl A 3 - coefrl (3) *rl A 2 - coef rl (4 ) *rl A l; 
atmp « atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coefr2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
atmp = atmp - coefr3 (1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coefr3 (4 ) *pxl A l; 
brl = [brl; atmp - coef r4 ( l)^pyl A 4 - coef r4 (2 ) *pyl A 3 - coef r4 ( 3) *pyl A 2 - coef r4 ( 4 ) *pyl A 
l;]; 

coefrS = pinv(arl) *real (brl) ; 

disp( 1 coefr5. . . 1 ) ; 

syms rl thl pxl pyl pzl; 

• llap = coefrl (1) *rl A 4 + coef rl (2) *rl A 3 + coef rl (3) *rl A 2 + coefrl (4) *rl A l; 
llap = rllap + coef r2 (1) *thl A 4 + coef r2 (2) *thl A 3 + coef r2 (3) *thl A 2 + coef r2 (4 ) *thl A l ; 
rllap = rllap + coef r3 (1) *pxl A 4 + coef r3 (2 ) *pxl A 3 + coef r3 (3) *pxl A 2 + coef r3 (4 ) *pxl A l; 
rllap = rllap + coefr4 (1) *pyl A 4 + coef r4 (2) *pyl A 3 + coefr4 (3) *pyl A 2 + coef r4 (4 ) *pyl A l; 
rllap = rllap- + coefr5 (1) *pzl A 4 + coef r5 (2) *pzl A 3 + coef r5 (3) *pzl A 2 + coefrS (4 ) *pzl A l; 
disp( 'rllap. . . 1 ) ; 

syms dBxn dByn dBzn dBxnap dBynap dBznap; 

syms diffy2 diffx2 diffz2 diffy2ap diffxap diffz2ap; 

diffy2 = evaKdiff (y2) ); 

diffx2 = diff (x2) ; 

diffz2 = evaKdiff (z2) ) ; 

xO = 0; 

yO = 0; 

zO = 0; 

rl = 1; 

pxl = 2.0; 

pyl = 0.0; 

pzl = 0.0; 

thl = 3.14159/2; 

arl = [thl A 4 thl A 3 thl A 2 thl A l;]; 
brl = [eval(diffy2) ;] ; 
thl = 3.14159; 

arl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 
brl = [brl; eval (dif f y2 ) ; ] ; 
thl = 3*3.14159/2; 

arl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 
^Drl = [brl; eval (dif fy2 );] ; 
^Bhl = 2*3.14159; 

£irl - [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 

brl = [brl; eval (dif fy2) ;] ; 

coefrl = pinv(arl) *brl; 

rl = 1; 



pxl = 0.0;. - 
pzl = 0.0; 
pzl = 0.0; 
thl = 0.0,- 

izl = [rl A 4 rl A 3 rl A 2 rl A l;]; 

5rl = [eval (diffy2) - coefrl (1) *thl A 4 - coef rl (2) *thl A 3 - coef rl (3) *thl A 2 - coefrl (4 ) *thl A 

1;]; 

rl - 2; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 

brl = [brl; eval (dif fy2) - cuef rl ( 1 ) *thl A 4 - coef rl (2) *thl A 3 - coefrl (3) *thl A 2 - coefrl (4) 
*thl A l; ] ; 
rl = 8; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 

brl = [brl; eval (dif fy2) - coef rl (!) *thl A 4 - coefrl (2) *thl A 3 - coefrl (3) *thl A 2 - coefrl (4) 

*thl A l;]; 

rl = 32; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 

brl = [brl; eval(diffy2) - coef rl ( 1) *thl A 4 - coef rl (2 ) *thl A 3 - coef rl (3) *thl A 2 - coefrl (4) 
*thl A l;]; 

coefr2 = pinv (arl) *brl; 

diffy2ap = coefrl (1) *thl A 4 + coef rl (2) *thl A 3 + coef rl (3) *thl A 2 + coefrl (4) *thl A l; 
diffy2ap = diffy2ap + coef r2 (1) *rl A 4 + coef r2 (2) *rl A 3 + coefr2 (3) *rl A 2 + coef r2 (4 ) *rl A l; 



xO = 


= 0; 




yO = 


= 0; 




z0 = 


= 0; 




rl = 


= 1; 




thl 


= 3.14159/2; 




arl 


= [thl A 4 thl A 3 thl A 2 thl A 


l;]; 


brl 


= [eval (dif fx2) ;] ; 




thl 


= 3.14159; 




arl 


= [arl; thl A 4 thl A 3 thl A 2 


thl A l;]; 


|rl 


= [brl; eval (dif f x2) ;] ; 




Phi 


« 3*3.14159/2; 




arl 


= [arl; thl A 4 thl A 3 thl A 2 


thl A l;]; 


brl 


= [brl; eval (dif f x2) ;] ; 




thl 


= 2*3.14159; 




arl 


= [arl; thl A 4 thl A 3 thl A 2 


thl A l;]; 


brl 


= [brl; eval (dif fx2) ;] ; 




coefrl = pinv(arl) *brl; 




rl = 


= 1; 




arl 


= [rl A 4 rl A 3 rl A 2 rl A l;]; 




brl 


= [eval (dif fx2) - coef rl ( 1 ) *thl A 4 



i;]; 
rl = 2; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 

brl « [brl; eval (dif fx2) - coefrl (1) *thl A 4 - coefrl (2) *thl A 3 - coef rl (3) *thl A 2 - coefrl (4) 
*thl A l;] ; 
rl = 8; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 

brl = [brl; eval(diffx2) - coef rl ( 1 ) *thl A 4 - coefrl (2) *thl A 3 - coefrl (3) *thl A 2 - coefrl (4) 

*thl A l;]; 

rl = 32; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 

brl = [brl; eval (dif fx2) - coefrl (1) *thl A 4 - coef rl (2) *thl A 3 - coef rl (3) *thl A 2 - coefrl (4) 
*thl A l;]; 

coefr2 = pinv (arl) *brl; 

diffx2ap = coefrl (1) *thl A 4 + coef rl (2 ) *thl A 3 + coef rl ( 3) *thl A 2 + coefrl (4 ) *thl A l; 
diffx2ap = diffx2ap + coef r2 (1) *rl A 4 + coefr2 (2) *rl A 3 + coef r2 (3) *rl A 2 + coefr2 (4) *rl A l; 
xO = 0; 

[0 - 0; 

\0 = 0; 
"rl = 1; 

thl = 3.14159/2; 

arl = [thl A 4 thl A 3 thl A 2 thl A l;]; 
brl = [diffz2;]; 



thl = 3.14159; 

arl = [arl; thl A 4thl A 3 thl A 2 thl A l;]; 
brl = [brl; diffz2;] ; 
thl = 3*3.14159/2; 

|rl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 
"rl = [brl; diffz2;] ; 
thl = 2*3.14159-; 

arl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 
brl = [brl; diffz2;] ; 
coefrl = pinv (arl) *brl; 
rl - 1; 

arl = [rl A 4 rl A 3 rl A 2 rl A l;]; 

brl = [diffz2 - coefrl (1) *thl A 4 - coefrl (2) *thl A 3 - coefrl (3) *thl A 2 - coefrl (4 ) *thl A l; ] ; 
rl = 2; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 

brl = [brl; diffz2 - coefrl (1) *thl A 4 - coefrl (2) *thl A 3 - coefrl (3) *thl A 2 - coefrl ( 4 ) *thl A 

l;]; 

rl = 8; 

arl - [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 

brl = [brl; diffz2 - coefrl (1) *thl A 4 - coefrl (2 ) *thl A 3 - coefrl ( 3) *thl A 2 - coefrl (4 ) *thl A 

i;]; 

rl = 32; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 

brl = [brl; diffz2 - coefrl (1) *thl. A 4 - coefrl (2) *thl A 3 - coefrl (3) *thl A 2 - coefrl (4 ) *thl A 
1;]; 

coefr2 = pinv (arl) *brl; 

diffz2ap = coefrl (l)*thl A 4 + coefrl (2) *thl A 3 + coefrl (3) *thl A 2 + coef rl ( 4 ) *thl A l; 
diffz2ap = diffz2ap + coef r2tD *rX"% + coef r2 (2) *rl A 3 + coefr2 (3) *rl A 2 + coef r2 (4 ) *rl A l; 



syms rl- thl pxl pyl pzl; 

dBxn = eval (dif f y2ap*rlz - rly*diffz2ap.) v ;/:%. i . . r 

~ ( Byn = eval (dif fx2ap*rlz - rlx*dif f z2ap) ;' 7 "'' ' 

Bzn = eval (dif fx2ap*rly rlx*dif fy2ap) ; 

%;arl = [rl pxl pyl pzl] : 

thl =0.0; 
pxl = 2.0; 
pyl = 0.0; 
pzl = 0.0; 
rl = 2; 

arl = [rl A 4 rl A 3 rl A 2 rl A l;]; 
brl = [eval (dBxn) ;] ; 
rl = 8; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 
brl = [brl; eval (dBxn) ;] ; 
rl = 16; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 
brl = [brl; eval (dBxn) ;] ; 
rl = 48; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 

brl = [brl; eval (dBxn) ;] ; 

coefrl = pinv (arl) *real (brl) ; 

disp ( 1 coefrl. . . 1 ) ; 

pxl = 0.0; 

pyl = 0.0; 

pzl = 0.0; 

rl = 1.0; 

thl = 3.14159/2; 

arl = [thl A 4 thl A 3 thl A 2 thl A l;]; 

jrl = [eval (dBxn) - coef rl ( 1 ) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; ] ; 
Phi = 3.14159; 
^arl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 

brl = [brl; eval (dBxn) - coef rl ( 1 ) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl ( 3) *rl A 2 - coef rl (4) *rr 

i;J; 

thl = 3*3.14159/2; 



arl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 

brl = [brl; eval(dBxn) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coefrl (3) *rl A 2 - coef rl (4 ) *rl A 
l;]; 

thl = 2*3.14159; 

•rl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 
rl = [brl; eval (dBxn) - coef rl ( 1) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4 ) *rl A 
1;]; 

coefr2 = pinv(arl) *real (brl) ; 

disp( 1 coefr2. . . 1 ) ; 

pyl = 0.0; 

pzl = 0.0; 

rl = 1.0; 

thl - 0.0; 

pxl = 2; 

arl = [pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = eval (dBxn) - coef rl ( 1 ) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
brl = [atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4) *thl A l; ] ; 
pxl = 8; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = eval (dBxn) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coefrl (3) *rl A 2 - coef rl (4 ) *rl A l; 
brl = [brl; atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4 ) *thl A 
l;]; 

pxl = 18; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = eval (dBxn) - coefrl (1) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
brl = [brl; atmp - coefr2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4) *thl A 
l;]; 

pxl = 4 8; t: 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = eval (dBxn) - coef rl ( 1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
brl = [brl; atmp - coefr2 (iy*-thl A 4 - coefr£(2) *thl A 3 - coef r2 ( 3) *thl A 2 - coef r2 (4 ) *thl A 

l;]; 

•oefr3 = pinv (arl) *real (brl) ; 
isp ( 1 coefr3. . . 1 ) ; 
pxl = 0.0; 
pzl = 0.0; 
thl = 0.0; 
rl = 0.0; 
pyl = 2; 

arl = [pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = eval (dBxn) - coefrl (1) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
atmp = atmp - coef r2 ( 1) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
brl = [atmp - coefr3 (1) *pxl A 4 - coef r3 (2) *pxl A 3 - coefr3 (3) *pxl A 2 - coefr3 (4) *pxl A l; ] ; 
pyl = 8; 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = eval (dBxn) - coef rl ( 1) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4) *thl A l; 
brl = [brl; atmp - coefr3 (1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 ( 4 ) *pxl A 
l;]; 

pyl =18; 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = eval (dBxn) - coef rl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coefr2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
brl = [brl; atmp - coefr3 (1) *pxl A 4 - coefr3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A 

i;]; 

pyl = 48; 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = eval (dBxn) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coefr2 (3) *thl A 2 - coefr2 (4) *thl A l; 
brl = [brl; atmp - coefr3 (1) *pxl A 4 - coefr3 (2) *pxl A 3 - coefr3 (3) *pxl A 2 - coef r3 (4 ) *pxl A 
^;]; 

flBoefr4 = pinv(arl) *real (brl) ; 
disp ( 1 coef r4 . . . 1 ) ; 
pxl = 0.0; 
pyl = 0.0; 
thl =0.0; 



rl = 0.0; 
pzl = 2; 

arl = [pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = eval(dBxn) - coefrl (1) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - coef rl ( 4 ) *rl A l ; 

•tmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coefr2 (3) *thl A 2 - coefr2 (4 ) *thl A l; 
tmp = atmp - coefr3 (1) *pxl A 4 - coefr3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
brl = [atmp - coefr4 (1) *pyl A 4 - coef r4 (2 ) *pyl A 3 - coefr4 (3) *pyl A 2 - coef r4 (4 ) *pyl A l; ] ; 
pzl = 8; " 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = eval(dBxn) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4) *thl A l; 
atmp = atmp - coefr3 (1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
brl = [brl; atmp - coef r4 ( 1) *pyl A 4 - coefr4 (2) *pyl A 3 - coef r4 (3) *pyl A 2 - coefr4 (4 ) *pyl A 

l;]; 

pzl = 32; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = eval(dBxn) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4) *thl A l; 
atmp = atmp - coefr3 (1) *pxl A 4 - coefr3 (2) *pxl A 3 - coefr3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
brl = [brl; atmp - coef r4 (1 ) *pyl A 4 - coef r4 (2 ) *pyl A 3 - coef r4 (3) *pyl A 2 - coef r4 (4 ) *pyl A 
1;]; 

pzl = 48; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp - eval(dBxn) - coef rl ( 1 ) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *th± A 4 - coef r2 (2) *th^L A 3- - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
atmp = atmp - coefr3 (1) *pxl A 4 - coef r3 (2^*pxl A .3 - coefr3 (3) *pxl A 2 - coef r3 ( 4 ) *pxl A l ; 
brl = [brl; atmp - coef r4 ( 1 ) *pyl A 4 - coefi:4 (2)l*pyl A 3 - coef r4 (3) *pyl A 2 - coef r4 (4 ) *pyl A 
l;]; | ..";> 

coefr5 = pinv(arl) *real (brl) ; } 
disp( f coefrS. . . 1 ) ; ' j ^ ■. 

syms rl thl pxl pyl pzl; **£"V\- 

dBxnap = coefrl (l)*rl A 4 + coefrl (2) *rl A 3- +>:Coef rl (3) * r l- A 2--+- coefrl (4 ) *rl A l; - 

•Bxnap = dBxnap + coefr2 (1) *thl A 4 + coef r2 ('2') *thl A 3 + coef r2 (3) *thl A 2 + coef r2 (4 ) *thl A I; 
Bxnap = dBxnap + coef r3 (1 ) *pxl A 4 + coefr3 (2) *pxl A 3 + coef r3 (3) *pxl A 2 + coefr3 (4 ) *pxl A l; 
dBxnap « dBxnap + coefr4 (1) *pyl A 4 + coefr4 (2) *pyl A 3 + coef r4 (3) *pyl A 2 + coef r4 (4 ) *pyl A l; 
dBxnap = dBxnap + coef r5 (1 ) *pzl A 4 + coef r5 (2) *pzl A 3 + coef r5 (3) *pzl A 2 + coef r5 (4 ) *pzl A l; 
disp ( 'dBxnap. 1 ) ; , 

syms rl thl pxl pyl pzl; 

rl = 1.0; 
pxl = 2.0; 
pyl =0.0; 
pzl = O.C- 
arl = [rl A 4 rl A 3 rl A 2 rl A l;]; 
brl = [eval (dByn) ;] ; 
rl = 8; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 
brl = [brl; eval (dByn) ;] ; 
rl = 16; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 
brl = [brl; eval (dByn) ;] ; 
rl = 48; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 
brl = [brl; eval (dByn) ;] ; 
coefrl = pinv(arl) *real (brl) ; 
disp( 'coefrl. . . ' ) ; 
thl = 3.14159/2; 

arl = [thl A 4 thl A 3 thl A 2 thl A l;]; 

brl = [eval (dByn) - coefrl (1) *rl A 4 - coef rl (2 ) *rl A 3 - coefrl (3) *rl A 2 - coef rl ( 4 ) *rl A l ; ] ; 
^^hl = 3.14159; 

fl|rl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 

^>rl = [brl; eval (dByn) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - coef rl ( 4 ) *rl A 
1;]; 

thl = 3*3.14159/2; 

arl o [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 



arl = 
fcrl = 



brl = [brl; eval(dByn) - coefrl (1) *rl A 4 - coef rl (2 ) *rl A 3 - coefrl (3) *rl A 2 - coefrl (4 ) *rl' 
i;]; 

thl = 2*3.14159; 

[arlr ttrf A 4 thl A 3 thl A 2 thl A l;]; 

[brl; eval(dByn) - coef rl ( 1 ) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4) *rl' 

coefr2 = pinv(arl) *real (brl) ; 
disp( f coefr2. . . 1 ) ; 
pyl =0.0 
pzl = 0.0 
thl = 0.0 
rl = 1.0; 
pxl = 2; 

arl = [pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = eval(dByn) - coefrl (1) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4) *rl A l; 
brl = [atmp - coefr2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 ( 4 ) *thl A l; ] ; 
pxl = 8; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 
atmp = eval(dByn) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - 
brl = [brl; atmp - coef r2 (1) *thl A 4 - coefr2 (2) *thl A 3 

l;]; 

pxl = 18; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 
atmp = eval(dByn) - coefrl (1) *rl A 4 - coef rl (2 ) *rl A 3 - 
brl = [brl; atmp - coef r2 (l)-*thl A 4 - coef r2 (2 ) *thl A 3 



coefrl (3) *rl A 2 - coefrl (4 ) *rl A l; 
coefr2(3)*thl A 2 - coef r2 (4 ) *thl' 



coefrl (3) *rl A 2 - coef rl {4 ) *rl A l; 
- coefr2 (3) *thl A 2 - coef r2 (4 ) *thl' 



pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 



coefrl (2) *rl A 3 - 
coefr2(2)*thl A 3 



coefrl (3) *rl A 2 - 
- coefr2 (3) *thl A 2 



coefrl (4) *rl A l; 
- coefr2 (4) *thl' 



l;]; 

pxl = 48; 
arl = [arl; t 

atmp = eval(dByn) - coef rl (1) *rl A 4 
brl = [brl; atmp - coefr2 (1) *thl A 4 

1;]; ~Z ^ 
coefr3 = pinv (arl ) *real (brl ) ; 
,isp( f coefr3. . . 1 ) ; 
Pxl = 0.0$ 
pzl = 0.0; 
thl = 0.0; 
rl = 1.0; 
pyl = 2; 

arl = [pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = eval(dByn) - coef rl ( 1 ) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coefr2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
brl = [atmp - coef r3 ( 1 ) *pxl A 4 - coefr3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; ] ; 
pyl = 8; 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = eval(dByn) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
brl = [brl; atmp - coefr3 (1) *pxl A 4 - coef r3 (2) *pxl A 3 - coefr3 (3) *pxl A 2 - coefr3 (4) *pxl A 
l;]; 

pyl = 18; 

arl - [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = eval(dByn) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 
atmp = atmp - coef r2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 



- coefrl (4) *rl A l; 
coefr2(4)*thl A l; 



brl = [brl; atmp - coefr3 (1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coefr3 (3) *pxl A 2 - coef r3 (4 ) *pxl 
l;]; 

pyl =48; 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = eval(dByn) - coef rl (1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4 
atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - 



rl A l; 
coefr2(4)*thl A l; 



brl = [brl; atmp - coef r3 (1 ) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3 ) *pxl A 2 - coefr3 (4 ) *pxl' 
l;]; 

joefr4 = pinv (arl) *real (brl) ; 
lisp( f coefr4 . . . ? ) ; 
pxl = 0.0; 
pyl = 0.0; • 
thl = 0.0; 
' rl = 1.0; 



pzl = 2; 

arl = [pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = eval(dByn) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 

•tmp = atmp - coefr3 (1) *pxl A 4 - coefr3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
rl = [atmp - coefr4 (1) *pyl A 4 - coefr4 (2) *pyl A 3 - coef r4 (3) *pyl A 2 - coef r4 ( 4 ) *pyl A l; ] ; 
pzl = 8; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = eval(dByn) - coef rl ( 1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
atmp = atmp - coef r3 (1) *pxl A 4 - coef r3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coefr3 (4) *pxl A l; 
brl = [brl; atmp - coef r4 (1) *pyl A 4 - coef r4 (2 ) *pyl A 3 - coef r4 (3) *pyl A 2 - coef r4 ( 4 ) *pyl A 
1;]; 

pzl = 32/ 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = eval(dByn) - coef rl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
atmp - atmp - coefr3 (1) *pxl A 4 - coefr3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
brl = [brl; atmp - coef r4 ( 1 ) *pyl A 4 - coef r4 (2 ) *pyl A 3 - coef r4 (3) *pyl A 2 - coef r4 (4 ) *pyl A 

i;]; 

pzl = 48; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = eval(dByn) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
atmp = atmp-- cuef r3 ( 1 ) *pxl~Z - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
brl = [brl; atmp - coef r4 ( 1 ) *pyl A 4 - coef r4 (2 ) *pyl A 3 - coefr4 ( 3) *pyl A 2 - coef r4 ( 4 ) *pyl A 
1;]; 

coefr5 = piny (arl) *real (brl fV 

disp( 1 coefrS. . . 1 ) ; 

syms rl thl pxl pyl pzl; 

dBynap = coef rl ( 1) *rl A 4 + coef rl (2 ) *rl A 3 + coef rl (3) *rl A 2 + coef rl (4 ) *rl A l; 

dBynap = dBynap + coef r2 ( 1) *thl A 4 + coef r2 (2) *thl A 3 + coef r2 (3) *thl A 2 + coef r2 (4 ) *thl A l 
j^Bynap = dBynap + coef r3 ( 1 ) *pxl A 4 + coef r3 (2) *pxl A 3 + coef r3 (3) *pxl A 2 + coef r3 (4 ) *pxl A l 
^PBynap = dBynap + coef r4 ( 1 ) *pyl A 4 + coef r4 ( 2 ) *pyl A 3 + coef r4 ( 3 ) *pyl A 2 + coef r4 ( 4 ) *pyl A l 

dBynap = dBynap + coef r5 (1) *pzl A 4 + coef r5 (2 ) *pzl A 3 + coef r5 (3) *pzl A 2 + coef r5 (4 ) *pzl A l 

disp ( 1 dBynap. . . 1 ) ; 



syms rl thl pxl pyl pzl; 



pxl = 2.0; 
pyl = 0.0; 
pzl = 0.0; 
thl = 0.0; 
rl = 1.0; 

arl = [rl A 4 rl A 3 rl A 2 rl A l;]; 
brl = [dBzn; ] ; 
rl = 8; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 
brl = [brl; dBzn; ] ; 
rl = 16; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 
brl = [brl; dBzn; ] ; 
rl = 48; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 

brl = [brl; dBzn;]; 

coefrl = pinv (arl) *real (brl) ; 

disp ( 1 coefrl . . . 1 ) ; 

pxl = 0.0; 

pyl = 0.0; 

pzl = 0.0; 

:1 = 1.0; 

|hl = 3.14159/2; 
^arl = [thi A 4 thl A 3 thl A 2 thl A l;]; 

brl = [dBzn - coef rl (1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4) *rl A l; ] ; 

thl = 3.14159; . 

arl = [arl; thl A 4 thl A 3 ^hl A 2 thl A l;]; 



brl = [brl; dBzn - coefrl (l)'*rl A 4 - coef rl (2) *rl A 3 
thl = 3*3.14159/2; 

arl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 
brl = [ferl-r dBzir - coef rl (If *rl A 4 - coefrl (2) *rl A 3 

•hi = 2*3~.14159; 
rl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 
brl = [brl; dBzn - coef rl (l)*rl A 4 - coef rl (2) *rl A 3 
coefr2 = pinv(arl) *real (brl) ; 
disp( 1 coefr2. . . 1 ) ; 
pyl = 0.0; 
pzl = 0.0; 
thl = 0.0; 
rl = 1.0; 
pxl = 2; 

arl = [pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 
atmp = dBzn - coef rl ( 1 ) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4 ) *rl A l; 
brl = [atmp - coefr2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 ( 3) *thl A 2 - coefr2 (4 ) *thl A l; ] ; 
pxl = 8; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = dBzn - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 

brl = [brl; atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4) *thl A 

l;]; 

pxl = 18; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = dBzn - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coefrl (3) *rl A 2 - coefrl (4) *rl A l; 

brl = [brl; atmp - coef r2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 ( 4 ) *thl A 

l;]; 

pxl = 48r; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = dBzn - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4) *rl A l; 

brl = [brl; atmp - coef r2 (1 ) : *thl A 4 - coef r2 (2 ) *thl A 3 - coefr2 (3) *thl A 2 - coefr2 (4) *thl A 

l;]; 

•oefr3 = pinv(arl) *real (brl) ; 
isp( 'coeft3. . . *~) ; 
pxl = 0.0; 
pzl = 0.0; 
thl = 0.0; 
rl = 1.0; 
pyl = 2; 

arl - [pyl A 4 pyl A 3 pyl A 2 pyT~l;]; 

atmp = dBzn - coef rl (1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
brl = [atmp - coef r3 (1 ) *pxl A 4 - coef r3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 ( 4 ) *pxl A l; ] ; 
pyl = 8; " 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = dBzn - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
brl = [brl; atmp - coefr3 (1) *pxl A 4 - coef r3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 ( 4 ) *pxl A 

l;]; 

pyl = 18; 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp « dBzn - coef rl (1 ) *rl A 4 - coef rl (2) *rl A 3 - coefrl (3) *rl A 2 - coef rl (4 ) *rl A l; 
atmp = atmp - coef r2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
brl = [brl; atmp - coef r3 (1) *pxl A 4 - coefr3 (2) *pxl A 3 - coefr3 (3) *pxl A 2 - coef r3 (4 ) *pxl A 

l;]; 

pyl = 48; 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = dBzn - coef rl (1) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4 ) *rl A l; 
atmp = atmp - coef r2 (1 ) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4 ) *thl A l; 
brl = [brl; atmp - coef r3 ( 1 ) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A 

^p oe fr4 = pinv (arl) *real (brl) ; 
disp( 1 coefr4 . . . 1 ) ; 
pxl = 0.0; 
pyl =0.0; 
rl = 1.0; 



- coefrl (3) *rl A 2 

- coefrl (3) *rl A 2 

- coefrl (3) *rl A 2 



coefrl (4) *rl A l;] ; 
coefrl (4) *rl A l;] ; 
coefrl (4) *rl A l;]; 



thl = 0.0; 
pzl = 2; 

arl = [pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp" « dBzn - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coefrl (3) *rl A 2 - coef rl (4 ) *rl A l; 

•tmp- = atmp - coef r2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coefr2 (3) *thl A 2 - coefr2 (4) *thl A l; 
tmp = atmp - coefr3 (1) *pxl A 4 - coefr3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
brl = [atmp - coefr4 (1) *pyl A 4 - coef r4 (2) *pyl A 3 - coefr4 (3) *pyl A 2 - coef r4 (4 ) *pyl A l; ] ; 
pzl = 8; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = dBzn - coef rl (1) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4) *rl A l; 
atmp = atmp - coef r2 ( 1 ) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4) *thl A l; 
atmp = atmp - coefr3 (1) *pxl A 4 - coefr3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coefr3 (4 ) *pxl A l; 
brl = [brl; atmp - coefr4 (1) *pyl A 4 - coef r4 (2) *pyl A 3 - coef r4 (3) *pyl A 2 - coef r4 (4 ) *pyl A 
1;]; 

pzl = 32; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = dBzn - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 - coefrl (4) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4 ) *thl A l; 
atmp = atmp - coefr3 (1) *pxl A 4 - coefr3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
brl = [brl; atmp - coef r4 ( 1) *pyl A 4 - coef r4 (2 ) *pyl A 3 - coefr4 (3) *pyl A 2 - coef r4 (4 ) *pyl A 
1;]; 

pzl = 48; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = dBzn - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - coef rl (4 ) *rl A l; 
atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coefr2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
atmp = atmp - coefr3 ( 1) *pxl A 4 - coef r3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
brl = [brl; atmp - coef r4 ( 1) *pyl A 4 - coef r4 (2) *pyl A 3 - coefr4 (3) *pyl A 2 - coef r4 ( 4 ) *pyl A 
1;]; 

coefrS = pinv(arl) *real (brl) ; 

disp( f coefr5. . . 1 ) ; 

syms rl thl pxl pyl pzl; 

dBznap = coefrl (1) *rl A 4 + coef rl (2) *rl A 3 + coef rl (3) *rl A 2 + coef rl (4 ) *rl A l; 
^dBznap = dBznap + coefr2 ( 1) *thl A 4 + coef r2 (2) *thl A 3 + coef r2 (3) *thl A 2 + coef r2 (4 ) *thl A l; 
^PBznap = dBznap + coef r3 ( 1 ) *pxl A 4 + coef r3 (2 ) *pxl A 3 + coef r3 (3) *pxl A 2 + coef r3 (4 ) *pxl A l; 

dBznap = dBznap + coefr4 (1) *pyl A 4 + coef r4 (2) *pyl A 3 + coef r4 (3) *pyl A 2 + coef r4 (4 ) *pyl A l; 

dBznap = dBznap + coef r5 (1) *pzl A 4 + coefr5 (2) *pzl A 3 + coef r5 (3) *pzl A 2 + coef r5 (4 ) *pzl A l ; 

disp ( 1 dBznap. . . 1 ) ; 

syms rl thl pxl pyl pzl; 

rl = 1.0; 
thl = 0.0; 
pxl = 2.0; 
pyl = 0.0; 
pzl =0.0; 

arl = [rl A 4 rl A 3 rl A 2 rl A l;]; 
brl = [real (eval (dBxnap) /eval (rllap) ) ; ] ; 
pxl = 8; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 
brl = [brl; real (eval (dBxnap) /eval (rllap) );] ; 
pxl = 16; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 
brl = [brl; real (eval (dBxnap) /eval (rllap) );] ; 
pxl = 48; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 

brl = [brl; real (eval (dBxnap) /eval (rllap) );] ; 

coefrl = pinv(arl) *real (brl) ; 

disp( 'coefrl. . . 1 ) ; 

pxl = 0.0; 

pyl = 0.0; 
^ozl = 0.0; 
■ l = 1.0; 
^thl = 3.14159/2; 

arl = [thl A 4 thl A 3 thl A 2 thl A l;]; 

brl = [real (eval (dBxnap) /e.val (rllap) ) - coefrl (1) *rl A 4 
coefrl (4) *rl A l;]; 



- coefrl (2) *rl A 3 - coef rl (3) *rl A 2 



thl = 3.14159; 

arl = [ari; thl A 4thl A 3 thl A 2 thl A l;]; 

brl = [brl; real (eval (dBxnap) /eval (rllap) ) - coef rl (1) *rl A 4 - coefrl (2) *rl A 3 - coefrl(3) 
*rl A 2 - coefrl (4 ) *rl A l; ] ; * 

•hi = 3*3.14159/2; 
rl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 
brl = [brl; real (eval (dBxnap) /eval (rllap) ) - coef rl ( 1) *rl A 4 - coef rl (2 ) *rl A 3 - coefrl (3) 
*rl A 2 - coefrl (4) *rl A l; ] ; 
thl = 2*3.14159; 

arl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 

brl = [brl; real (eval (dBxnap) /eval (rllap) ) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coefrl (3) 

*rl A 2 - coefrl (4) *rl A l; ] ; 

coefr2 = pinv (arl) *real (brl) ; 

disp ( 1 coefr2 ...'); 

pxl = 2.0; 

pyl - 0.0; 

thl = 0.0; 

rl = 1.0; 

arl = [pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = real (eval (dBxnap) /eval (rllap) ) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coefrl (3) *rl A 2 - 
coefrl (4) *rl A l; 

brl = [atmp - coefr2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4 ) *thl A l; ] ; 
pxl = 8; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = real (eval (dBxnap) /eval (rllap) ) - coef rl ( 1 ) *rl A 4 - coefrl (2) *rl A 3 - coefrl (3) *rl A 2 - 
coefrl (4) *rl A l; 

brl = [brl; atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 ( 4 ) *thl A 
i;]; .:T 

pxl = 18; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = rep (eval (dBxnap) /ev*C (rllap)) - coef rl ( 1 ) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - 
coefrl (4) *rl A l; 

^brl = [brl; atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A 
pxl = 48; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = real (eval (dBxnap) /eval (rllap) ) - coef rl ( 1 ) *rl A 4 - coef rl (2 ) *rl A 3 - coefrl (3) *rl A 2 - 
coefrl (4) *rl A l; 

brl = [brl; atmp - coef r2 (1 ) *thl A 4 - coef r2 (2) *thl A 3 - coefr2 (3) *thl A 2 - coefr2 (4 ) *thl A 
l;]; 

coefr3 = pinv (arl) *real (brl ) ; 

disp( f coefr3. . . 1 ) ; 

pxl = 0.0; 

pzl = 0.0; 

thl = 0.0; 

rl = 1.0; 

pyl = 2; 

arl = [pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp == real (eval (dBxnap) /eval (rllap) ) - coef rl ( 1 ) *rl A 4 - coef rl (2 ) *rl A 3 - coefrl (3) *rl A 2 - 
coefrl (4) *rl A l; 

atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
brl = [atmp - coefr3 (1) *pxl A 4 - coef r3 (2) *pxl A 3 - coefr3 (3) *pxl A 2 - coefr3 (4 ) *pxl A l; ] ; 
pyl = 8; 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp - real (eval (dBxnap) /eval (rllap) ) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coefrl (3) *rl A 2 - 
coefrl (4) *rl A l; 

atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l ; 
brl = [brl; atmp - coefr3 (1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 ( 4 ) *pxl A 

l;]; 

pyl = 18; 

^^rl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

^Btmp = real (eval (dBxnap) /eval (rllap) ) - coef rl ( 1 ) *rl A 4 - coef rl ( 2 ) *rl A 3 - coef rl ( 3) *rl A 2 - 

coefrl (4) *rl A l; 

atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
brl = [brl; atmp - coe_f x3 (,1 ) *pxl"4 - coef r3 (2 ) *pxl"3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A 
l;]; 



pyl = 48; 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = real(eval(dBxnap)/eval(rllap) ) - coef rl ( 1 ) *rl A 4 - coef rl (2 ) *rl A 3 - coefrl (3) *rl A 2 - 
coefrl (4) *rl A l; 

•tmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4 ) *thl A l; 
rl = [brl; atmp - coef r3 ( 1 ) *pxl A 4 - coef r3 (2) *pxl A 3 - coefr3 (3) *pxl A 2 - coef r3 (4 ) *pxl A 

coefr4 = pinv (arl) *real (brl) ; 

disp ( 1 coef r4 . . . 1 ) ; 

pxl = 0.0; 

pyl = 0.0; 

thl = 0.0; 

rl = 1.0;., 

pzl = 2; 

arl = [pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = real(eval(dBxnap)/eval(rllap) ) - coefrl (1) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - 
coefrl (4) *rl A l; 

atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4 ) *thl A l; 
atmp = atmp - coefr3 (1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
brl = [atmp - coefr4 ( 1) *pyl A 4 - coef r4 (2 ) *pyl A 3 - coef r4 (3) *pyl A 2 - coefr4 (4 ) *pyl A l; ] ; 
pzl = 8; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = real(eval(dBxnap) /eval(rllap) ) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - 
coefrl (4) *rl A l; 

atmp = atmp - coefr2 (1) *thl^4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4 ) *thl A l; 
atmp = atmp - coefr3 (1) *pxl A 4 - coef r3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coefr3 (4) *pxl A l; 
brl = [brl; atmp - coef r4 (1) *pyl A 4 - coefr4 (2) *pyl A 3 - coef r4 (3) *pyl A 2 - coefr4 (4) *pyl A 
1;]; 

pzl = 32; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = real(eval(dBxnap) /ev#l (rllarr) ) - coefrl (1) *rl A 4 - coef rl (2 ) *rl A 3 - coefrl (3) *rl A 2 - 
coefrl (4) *rl A l; 

•tmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coefr2 (3) *thl A 2 - coef r2 ( 4 ) *thl A l ; 
tmp = atmp - coefr3 (1 ) *pxl A 4 - coefr3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 ( 4 ) *pxl A l; 
brl = [brl; atmp - coefr4 (1) *pyl A 4 - coef r4 (2 ) *pyl A 3 - coef r4 (3) *pyl A 2 - coef r4 (4 ) *pyl A 

i;]; 

pzl = 48; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = real(eval(dBxnap) /eval(rllap) ) - coef rl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - 
coefrl (4) *rl A l; 

atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coefr2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
atmp = atmp - coefr3 (1) *pxl A 4 - coefr3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coefr3 (4) *pxl A l; 
brl = [brl; atmp - coefr4 (1) *pyl A 4 - coef r4 (2) *pyl A 3 - coef r4 (3) *pyl A 2 - coefr4 (4) *pyl A 

l;]; 

coefr5 = pinv(arl) *real (brl) ; 

disp ( 1 coef r5 . . . ? ) ; 

syms rl thl pxl pyl pzl; 

dBx = coefrl (1) *rl A 4 + coef rl (2 ) *rl A 3 + coef rl (3) *rl A 2 + coefrl (4 ) *rl A l; 
dBx = dBx + coefr2(l)*thl A 4 + coef r2 (2) *thl A 3 + coef r2 (3) *thl A 2 + coefr2 (4) *thl A l; 
dBx = dBx + coefr3(l)*pxl A 4 + coef r3 (2) *pxl A 3 + coef r3 (3) *pxl A 2 + coefr3 ( 4 ) *pxl A l; 
dBx = dBx + coefr4 (1) *pyl A 4 + coef r4 (2 ) *pyl A 3 + coef r4 (3) *pyl A 2 + coef r4 (4 ) *pyl A l; 
dBx = dBx + coefr5 (1) *pzl A 4 + coef r5 (2) *pzl A 3 + coef r5 (3) *pzl A 2 + coefr5 (4 ) *pzl A l; 
disp ( f dBx. . . 1 ) ; 



rl = 1.0; 

circlelx = u0*int (dBx, thl, 0, 2*3. 14159) ; 
disp ( ' finished circlelx. . . 1 ) ; 

^^printf (fid, 1 \n\nCirclelx = Sum of: \n f ); 
■or il = 1:4 

WV fprintf (fid, f %20. 12f» , coefrl (il) ) ; 
fprintf (fid, f pxl A %li' , il) ; 
fprintf (fid, ! \n' ) ; 
end 



for il = 1:4 Z 

fprintf (fid, 1 %20 . 12f f , coef r2 (il ) ); 
fprintf (fid, ' pyl A %li ' , il ) ; 

f print ftffdi h \n' ); " 

f d 

>r il = 1:4 

fprintf (fid, f %20 . 12f 1 , coef r3 (il ) ); 
fprintf (fid, ' pzl A %li ' , il ) ; 
fprintf (fid, f \n f ) ; 
end 

for il = 1:4 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, '\n' ) ; 
end 

for il = 1:4 

fprintf (fid, ' %20 . 12f 1 , coef r5 ( il ) ) ; 

fprintf (fid, 1 th2 A %li ' , il) ; 

fprintf (fid, '\n f ) ; 
end 

syms tcirclelx tcirclely tcirclelz testl test2 test3 test4 testS; 
tcirclelx = coefrl (1) *testl A 4 + coef rl (2) *testl A 3 + coef rl (3) *testl A 2 

+ coefr2 (1) *test2 A 4 



! %20.12f f ,coefr4 (il) ) ; 
• thl A %li f ,il) ; 



m 



tcirclelx = tcirclelx 
(4) *test2 A l; 
tcirclelx = tcirclelx 
(4) *test3 A l; 
tcirclelx = tcirclelx 
(4)*test4 A l; 
tcirclelx = tcirclelx 
(4) *test5 A l; 
testl = 2*rl; /%; * 
test 2 = 0.0; 
k est3 = 0.0; 
st4 = 0.0; 
tests = 0.0; 

fprintf (fid, 1 \n\nTest of 
fprintf (fid, ' \npxl = %20 
fprintf (fid, 1 \npyl = %20 
fprintf (fid, '\npzl = %20 
fprintf (fid, 1 \nthl = %20 
fprintf (fid, 1 \nth2 = %20 
fprintf (fid, 1 \ncirclelx = 
fprintf (fid, 'Xn 1 ) ; 



coefr3(l) *test3 A 4 
coefr4 (1) *test4 A 4 
coefr5 (1) *test5 A 4 



coefr2 (2) *test2 A 3 
coefr3(2) *test3 A 3 
coefr4 (2) *test4 A 3 
coefrS (2) *test5 A 3 



+ coefrl (4) *testl A l; 
coefr2 (3) *test2 A 2 + coefr2 



coefr3 (3) *test3 A 2 + coefr3 
coefr4 (3) *test4 A 2 + coefr4 
coefr5 (3) *test5 A 2 + coefr5 



circlelx: 1 ) 
12f 1 , testl) 
12f 1 , test2) 
12f ' , test3) 
12f 1 , test4) 
12f \ testS) 
%20.12f • ,eval (tcirclelx) ) ; 



syms rl thl pxl pyl pzl; 



rl = 


= 1.0; 




thl 


- 0.0; 




pxl 


- 0.0; 




pyl 


= 0.0; 




pzl 


= 0.0; 




arl 


= [rl A 4 


rl A 3 rl A 2 rl A l;] ; 


brl 


= [real (eval (dBynap) /eval (rllap) ) ; ] ; 


rl = 


= 8; 




arl 


= [arl; 


rl A 4 rl A 3 rl A 2 rl A l;] ; 


brl 


- [brl; 


real (eval (dBynap) /eval (rllap) ) ; ] ; 


rl = 


= 16; 




arl 


= [arl; 


rl A 4 rl A 3 rl A 2 rl A l;] ; 



brl = [brl; real (eval (dBynap) /eval (rllap) );] ; 
ff l = 48; 

Irl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 

brl = [brl; real (eval (dBynap) /eval (rllap) );] ; 

coefrl = pinv(arl) *real (brl) ; 

disp( 'coefrl. . . 1 ) ; 

pxl = 0.0; 



pyl = 0.0; 
thl = 0.0; 
rl = 1.0; 
thl = 20; 

•rl = [thl A 4 thl A 3 thl A 2 thl A l;]; 
rl = [real(eval(dBynap) /eval(rllap) ) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 
coefrl (4) *rl A l;]; 
thl = 80; 

arl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 

brl = [brl; real (eval (dBynap) /eval (rllap) ) - coefrl (1) *rl A 4 - coef rl (2 ) *rl A 3 - coefrl (3) 
*rl A 2 - coefrl (4) *rl A l;]; 
thl = 320; 

arl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 

brl = [brl; real (eval (dBynap) /eval (rllap) ) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coefrl (3) 
*rl A 2 - coefrl (4) *rl A l; ] ; 
thl = 720; 

arl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 

brl = [brl; real (eval (dBynap) /eval (rllap) ) - coef rl ( 1 ) *rl A 4 - coef rl (2) *rl A 3 - coefrl (3) 

*rl A 2 - coefrl (4) *rl A l;]; 

coefr2 = pinv (arl) *real (brl) ; 

disp( 1 coefr2. . . 1 ) ; 

pyl = 0.0; 

pxl = 0.0; 

rl = 1.0; 

thl = 0.0; 

pzl = 2; 

arl = [pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = real (eval (dBynap) /eval (rllap) ) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 
coefrl (4) *rl A l; 

brl - [atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; ] ; 
pzl = 8; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 
^^tmp = real (eval (dBynap) /eval (rllap) ) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 
^Poefrl (4) *rl A l; 

brl = [brl; atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4 ) *thl A 

l;]; 

pzl = 18; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = real (eval (dBynap) /eval (rllap) ) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 
coefrl (4) *rl A l; 

brl = [brl; atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A 
l;]; 

pzl = 48; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = real (eval (dBynap) /eval (rllap) ) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 
coefrl (4) *rl A l; 

brl = [brl; atmp - coef r2 ( 1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A 
l;]; 

coefr3 = pinv(arl) *real (brl) ; 

disp ( 1 coefr3. . . 1 ) ; 

pxl = 0.0; 

pzl = 0.0; 

rl = 1.0; 

thl = 0.0; 

pyl = 2; 

arl = [pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = real (eval (dBynap) /eval (rllap) ) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coefrl (3) *rl A 2 
coefrl (4) *rl A l; 

atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coefr2 (3) *thl A 2 - coefr2 (4) *thl A l; 
brl = [atmp - coefr3 (1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coefr3 (4) *pxl A l; ] ; 
^?yl = 8; 

^■rl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = real (eval (dBynap) /eval (rllap) ) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 
coefrl (4) *rl A l; 

atmp = atmp - coefr.2 (l)..tthl A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 ( 4 ) *thl A l; 
brl = [brl; atmp - coefr3 (1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coefr3 (4 ) *pxl A 



l;]; 

pyl = 18/ 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = real (evafl (dBynap) /eval (rllap) ) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 

•oefrl (4) *rl A l; 
tmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4 ) *thl A l; 
brl = [brl; atmp - coef r3 (1 ) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 ( 4 ) *pxl A 
l;]; 

pyl = 48; 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = real (eval (dBynap) /eval (rllap) ) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coefrl (3) *rl A 2 
coefrl (4) *rl A l; 

atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4) *thl A l; 
brl = [brl; atmp - coef r3 (1) *pxl A 4 - coef r3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coefr3 (4) *pxl A 
l;]; 

coefr4 = pinv(arl) *real (brl) ; 

disp ( 1 coefr4 . . . 1 ) ; 

pxl = 0.0; 

pyl = 0.0; 

rl = 1.0; 

thl = 0.0; 

pzl = 2; 

arl = [pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = real (eval (dBynap) /eval (rllap) ) - coef rl ( 1 ) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 
coefrl (4) *rl A l; 

atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
atmp = atmp - coef r3 ( 1 ) *pxl A 4 - coef r3 (2) *pxl A 3 - coef r3 ( 3) *pxl A 2 - coefr3 (4 ) *pxl A l; 
brl = [atjnp* - ccrefr4 (1) *pyl~4 - coe-fr4 (2) *pyl A 3 - coef r4 (3) *pyl A 2 - coef r4 ( 4 ) *pyl A l; ] ; 
pzl = 8; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = real (eval (dBynap) /evil (rllap) ) - coefrl ( 1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 
coefrl (4) *rl A l; 

^^tmp = atmp - coef r2 ( 1 ) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
^Ptmp = atmF - caefr3(l)*pxr*4 - coefr3 (2) *pxl A 3 - coefr3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
- brl = [brl; atmp - coefr4 (1) *pyl A 4 - coef r4 (2 ) *pyl A 3 - coefr4 (3) *pyl A 2 - coef r4 (4 ) *pyl A 
1;]; 

pzl = 32; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = real (eval (dBynap) /eval (rllap) ) - coef rl ( 1 ) *rl A 4 - coef rl (2 ) *rl A 3 - coefrl (3) *rl A 2 
coefrl (4 ) Jt1 a 1; 

atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coefr2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
atmp = atmp - coefr3 ( 1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 ( 4 ) *pxl A l; 
brl = [brl; atmp - coefr4 (1) *pyl A 4 - coef r4 (2 ) *pyl A 3 - coef r4 ( 3) *pyl A 2 - coef r4 (4 ) *pyl A 

l;]; 

pzl = 48; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = real (eval (dBynap) /eval (rllap) ) - coefrl (1) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 
coefrl (4) *rl A l; 

atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
atmp = atmp - coefr3 (1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coefr3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
brl = [brl; atmp - coef r4 (1) *pyl A 4 - coef r4 (2 ) *pyl A 3 - coefr4 (3) *pyl A 2 - coef r4 (4 ) *pyl A 

l;]; 

coefr5 = pinv (arl) *real (brl) ; 

disp( f coefr5. . . 1 ) ; 

syms rl thl pxl pyl pzl; 



dBy 


= coefrl (1) *rl A 4 + coefrl 


(2)*rl A 3 + coefrl ( 


3)*rl A 2 + 


coefrl (4) 


*rl A l; 




dBy 


= dBy 


+ 


coefr2(l)*thl A 4 


+ 


coefr2(2)*thl A 3 


+ 


coefr2 (3) 


*thl A 2 


+ 


coefr2 (4) 


*thl A l 


dBy 


= dBy 




coefr3 (1) *pxl A 4 


+ 


coefr3(2)*pxl A 3 


+ 


coefr3 (3) 


*pxl A 2 


+ 


coefr3 (4 ) 


*pxl A l 


dBy 


= dBy 


+ 


coefr4 (1) *pyl A 4 


+ 


coefr4 (2) *pyl A 3 


+ 


coefr4 (3) 


*pyl A 2 


+ 


coefr4 (4) *pyl A l 


dBy 


= dBy 


+ 


coefrS (1) *pzl A 4 


+ 


coefr5(2)*pzl A 3 


+ 


coefr5(3)*pzl A 2 


+ 


coefrS (4) 


*pzl A l 



^iisp ( 'dBy. . . 1 ) ; 
rl = 1.0; 



circlely = u0*int (dBy^.thl, 0, 2*3. 14159) ; 
disp ( 1 finished circl-ely. . . 1 ) ; 



fprintf (fid, 1 \n\nCirclely = Sum of: \n f ); 
for il = 1:4 

fprintf ( fid, ' %20 . 12f ' , coef rl (il ) ) ; 

fprintf (fid, 1 pxl A %li\ il)-; 

fprintf (fid, '\n' ) ; 
end 

for il = ,1:4 

fprintf (fid, 1 %20.12f ' , coefr2 (il) ) ; 

fprintf (fid, 1 pyl A %li 1 , il) ; 

fprintf (fid, ! \n' ) ; 
end 

for il = 1:4 

fprintf (fid, ' %20. 12f 1 , coefr3 (il) ) ; 

fprintf (fid, 1 pzl A %li f , il) ; 

fprintf (fid, f \n f ) ; 
end 

for il = 1:4 

fprintf (fid, 1 %20. 12f f , coefr4 (il) ) ; 

fprintf (fid, 1 thl A %li 1 , il) ; 

fprintf (fid, '\n' ) ; 
end 

for il = 1:4 

fprintf (fid, ' %20.12f f , coefr5 (il) ) ; 

fprintf (fid, 1 th2 A %li f , il) ; 

fprintf (fid, '\n' ) ; 
end 

syms testl test2 test3 test4 test5; 

tcirclely = coefrl ( 1 ) *testl A 4 + coef rl (2) *testl A 3 + coef rl (3) *testl A 2 + coefrl (4) *testl A l; 
tcirclely = tcirclely + coef r2 ( 1 ) *test2 A 4 + coef r2 (2) *test2 A 3 + coef r2 (3) *test2 A 2 + coefr2 
(4)*test2 A l; 

tcirclely = tcirclely + coef r3 ( 1 ) *test3 A 4 + coef r3 (2 ) *test3 A 3 + coef r3 (3 ) *test3 A 2 + coefr3 
^^4) *test3 A l; 

^^circlely = tcirclely + coef r4 ( 1 ) *test4 A 4 + coef r4 (2 ) *test4 A 3 + coef r4 ( 3 ) *test4 A 2 + coefr4 
(4)*test4 A l; 

tcirclely = tcirclely + coef r5 ( 1 ) *test5 A 4 + coefrS (2) *test5 A 3 + coefr5 (3) *test5 A 2 + coefr5 

(4) *test5 A l; 

testl - 2*rl; 

test2 = 0.0; 

test3 = 0.0; 

test4 = 0.0; 

test5 = 0.0; 

fprintf (fid, 1 \n\nTest of circlely: 1 ) ; 

fprintf (fid, 1 \npxl = %20 . 12f 1 , testl ) ; 

fprintf (fid, f \npyl = %20 . 12f f , test2 ) ; 

fprintf (fid, f \npzl = %20. 12f 1 , test3) ; 

fprintf (fid, 1 \nthl = %20 . 12f 1 , test4 ) ; 

fprintf (fid, 1 \nth2 = %20 . 12f f , test5) ; 

fprintf (fid, ! \ncirclely = %20 . 12f 1 , eval (tcirclely) ) ; 

fprintf (fid, '\n' ) ; 



syms rl thl pxl pyl pzl; 



pxl =0.0 
pyl =0.0 
pzl = 0.0 
thl = 0.0 
rl = 2; 

arl = [rl A 4 rl A 3 rl A 2 rl A l;]; 
Drl = [real (eval (dBznap) /eval (rllap) );] ; 
Jl = 8; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 

brl = [brl; real (eval (dBznap) /eval (rllap) );] ; 

rl = 16; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 



brl = [brl; rea.l (eval (dBznapr) /eval (rllap) );] ; 
rl = 48; 

arl = [arl; rl A 4 rl A 3 rl A 2 rl A l;]; 
brl = [brl; real (eval (dBznap) /eval (rllap) );] ; 
"^oefrl = pinv(arl) *real (brl) ; 
!isp( 'coefrl. . . 1 ) ; 
pxl = 0.0; 
pyl = 0.0; 
pzl = 0.0; 
rl = 1.0; 
thl = 3.14159/2; 

arl - [thl A 4 thl A 3 thl A 2 thl A l;]; 

brl = [real (eval (dBznap) /eval (rllap) ) - coefrl (1) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - 
coefrl(4)*rl A l;] ; 
thl = 3.14159; 

arl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 

brl = [brl; real (eval (dBznap) /eval (rllap) ) - coef rl (1) *rl A 4 - coef rl (2 ) *rl A 3 - coefrl(3) 
*rl A 2 - coefrl (4) *rl A l;] ; 
thl = 3*3.14159/2; 

arl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 

brl = [brl; real (eval (dBznap) /eval (rllap) ) - coef rl ( 1 ) *rl A 4 - coef rl (2) *rl A 3 - coefrl(3) 
*rl A 2 - coefrl (4) *rl A l; ] ; 
thl = 2*3.14159; 

arl = [arl; thl A 4 thl A 3 thl A 2 thl A l;]; 

brl = [brl; real (eval (dBznap) /eval (rllap) ) - coef rl (1 ) *rl A 4 - coef rl (2) *rl A 3 - coefrl(3) 

*rl A 2 - coefrl(4)*rl A l;] ; 

coefr2 = pinv (arl ) *real (brl ) ; 

disp ( 1 coefr2 ...'); 

pyl = 0.0; 

pzl = 0.0; 

rl = l.Or. 

thl = 0.0; 

axl * 2; 

^|rl = [pxl A 4 pxl A 3 pxl A 2 prf A l;]; 

atmp = real (eval (dBznap) /eval (rllap) ) - coefrl (1) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - 
coefrl(4)*rl A l; 

brl = [atmp - coefr2 (1) *thl*4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; ] ; 
pxl = 8; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = r earl (eval (dBznap) /evstl(rl lap) ) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coefrl (3) *rl A 2 - 
coefrl (4) *rl A l; 

brl = [brl; atmp - coefr2 (1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 ( 4 ) *thl A 

l;]; 

pxl = 18; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp = real (eval (dBznap) /eval (rllap) ) - coef rl ( 1 ) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 - 
coefrl (4) *rl A l; 

brl = [brl; atmp - coefr2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coefr2 (3) *thl A 2 - coefr2 (4) *thl A 

l;]; 

pxl = 48; 

arl = [arl; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

atmp « real (eval (dBznap) /eval (rllap) ) - coef rl ( 1) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 - 
coefrl (4) *rl A l; 

brl = [brl; atmp - coefr2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A 

l;]; 

coefr3 = pinv(arl) *real (brl) ; 

disp( 1 coefr3. . . ' ) ; 

pxl = 0.0; 

pzl = 0.0; 

rl = 1.0; 

-hi = 0.0; 

)yl = 2; 

"arl = [pyl A 4 pyl A 3 pyl^2 pyl A l;]; 

atmp = real (eval (dBznap) /eval (rllap) ) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coefrl (3) *rl A 2 - 
coefrl (4) *rl A l; 

atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 



brl = [atmp - coefr3 (1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; ] ; 
pyl = 8; 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = real (eval (dBznap) /eval (rllap) ) - coef rl ( 1 ) *rl A 4 - coef rl (2 ) *rl A 3 - coefrl (3) *rl A 2 

•oefrl(4)*rl A l; 
tmp = atmp - coef r2 ( 1) *thl A 4 - coef r2 (2 ) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
brl = [brl; atmp - coef r3 (1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 ( 4 ) *pxl A 
l;]; 

pyl = 16; 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = real (eval (dBznap) /eval (rllap) ) - coefrl (1) *rl A 4 - coef rl (2) *rl A 3 - coef rl (3) *rl A 2 
coefrl (4) *rl A l; 

atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coefr2 (4 ) *thl A l; 
brl = [brl; atmp - coefr3 (1) *pxl A 4 - coefr3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coefr3 (4) *pxl A 
l;]; 

pyl = 64; 

arl = [arl; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

atmp = real (eval (dBznap) /eval (rllap) ) - coefrl (1) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 
coefrl (4) *rl A l; 

atmp = atmp - coefr2 (1) *thl A 4 - coef r2 (2) *thl A 3 - coefr2 (3) *thl A 2 - coefr2 (4) *thl A l; 
brl = [brl; atmp - coef r3 (1) *pxl A 4 - coef r3 (2) *pxl A 3 - coefr3 (3) *pxl A 2 - coef r3 (4 ) *pxl A 
l;]; 

coefr4 = pinv (arl) *real (brl) ; 

disp ( 1 coefr4 . . . 1 ) ; 

pxl = 0.0; 

pyl = 0.0; 

rl = 1.0; 

thl = 0.0.2 

pzl = 2; 

arl = [pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = real (eval (dBznap) /evarl (rllap) ) - coef rl (1) *rl A 4 - coefrl (2) *rl A 3 - coefrl (3) *rl A 2 
coefrl (4) *rl A l; 

^^tmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coefr2 (3) *thl A 2 - coefr2 (4 ) *thl A l; 
^Ptmp = atmp - coef r3 ( 1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coefr3 (4 ) *pxl A l; 

brl = [atmp - coef r4 ( 1) *pyl A 4 - coef r4 (2 ) *pyl A 3 - coef r4 (3) *pyl A 2 - coefr4 (4 ) *pyl A l; ] ; 

pzl = 8; 

arl = [arlr pzf^ pzl^J pzT*2 pzl A lr]; 

atmp = real (eval (dBznap) /eval (rllap) ) - coef rl ( 1 ) *rl A 4 - coefrl (2) *rl A 3 - coef rl (3) *rl A 2 
coefrl (4) *rl A l; 

atmp = atmp - coefr2 (1) *thi A 4 - coef r2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
atmp = atmp - coefr3 (1) *pxl A 4 - coef r3 (2) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
brl = [brl; atmp - coefr4 (1) *pyl A 4 - coef r4 (2 ) *pyl A 3 - coef r4 (3) *pyl A 2 - coef r4 (4 ) *pyl A 

l;]; 

pzl = 32; 

arl - [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp = real (eval (dBznap) /eval (rllap) ) - coefrl (1) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 
coefrl (4) *rl A l; 

atmp = atmp - coefr2 (1) *thl A 4 - coefr2 (2) *thl A 3 - coef r2 (3) *thl A 2 - coef r2 (4 ) *thl A l; 
atmp = atmp - coefr3 (1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coefr3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
brl « [brl; atmp - coefr4 (1) *pyl A 4 - coef r4 (2) *pyl A 3 - coef r4 (3) *pyl A 2 - coef r4 (4 ) *pyl A 
l;]; 

pzl = 48; 

arl = [arl; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

atmp - real (eval (dBznap) /eval (rllap) ) - coef rl ( 1 ) *rl A 4 - coef rl (2 ) *rl A 3 - coef rl (3) *rl A 2 
coefrl (4) *rl A l; 

atmp = atmp - coef r2 ( 1) *thl A 4 - coef r2 (2) *thl A 3 - coef r2 ( 3) *thl A 2 - coefr2 (4) *thl A l; 
atmp = atmp - coefr3 (1) *pxl A 4 - coef r3 (2 ) *pxl A 3 - coef r3 (3) *pxl A 2 - coef r3 (4 ) *pxl A l; 
brl = [brl; atmp - coef r4 ( 1 ) *pyl A 4 - coefr4 (2) *pyl A 3 - coef r4 ( 3) *pyl A 2 - coef r4 (4 ) *pyl A 
l;] ; 

coefr5 = pinv (arl) *real (brl ) ; 
^^isp ( 1 coef r5 ...'); 
^Byms rl thl pxl pyl pzl; 

ciBz = coefrl (1) *rl A 4 + coef rl (2 ) *rl A 3 + coefrl (3) *rl A 2 + coefrl (4 ) *rl A l; 

dBz - dBz,+ coefr2(l) *thl A 4 + coefr2 (2) *thl A 3 + coef r2 (3) *thl A 2 + coef r2 (4 ) *thl A l; 

dBz = dBz" + coef rltUJjRXl£A. + coefr3 (2) *pxl A 3 + coefr3 (3) *pxl A 2 + coef r3 (4 ) *pxl A l; 

dBz = dBz + coefr4^^h*S^ ; ^ + coef r4 (2 ) *pyl A 3 '+ coef r4 (3);*pyd A 2 + coef r4 (4 ) *pyl A l; 



dBz = dBz + coefr5(l) *pzl A 4 .J- coef r5 (2) *pzl A 3 + coef r5 (3) *pzl A 2 + coef r5 ( 4 ) *pzl A l; 
disp ( 'dBz'. . . 1 ) ; 



1 = l.Or 

circlelz = uO*int (dBz, thl, 0, 2*3. 14159) ; 
disp ( 1 finished circlelz . . . 1 ) ; 

fprintf (f id, 1 \n\nCirclelz = Sum of: \n'); 
for il = 1:4 

fprintf (fid, '%20.12f * , coefrl (il) ) ; 

fprintf (fid, 1 pxl A %li ' , il) ; 

f print f (fid, ' \n' ) ; 
end 

for il = 1:4 

f printf ( fid, ? %20 . 12f ' , coef r2 ( il ) ) ; 

fprintf (fid, ' pyl A %li f , il) ; 

f printf (fid, '\n' ) ; 
end 

for il = 1:4 

fprintf (fid, * %20. 12f 1 , coefr3 (il) ) ; 

fprintf (fid, 1 pzl A %li ' , il) ; 

fprintf (fid, ' \n' ) ; 
end 

for il = 1:4 

fprintf (fid, 1 %20 . 12f 1 , coef r4 (il) ) ; 

fprintf (fid, 1 thl A %li 1 , il) ; 

fprintf (fid, ' \n f ) ; 
end 

for il = 1:4 

f print ftfxd, -*%2 0 . 1 2 f^ , coefx5 ( i 1 ) ) ; 

fprintf (fid, 1 th2 A %li 1 , il) ; 

fprintf (fid, *\n') ; 
end ^ 

syms testl test2 test3 test4 test5; 

tcirclelz = coefrl (1) *testl A 4 + coef rl (2 ) *testl A 3 + coef rl (3) *testl A 2 + coefrl (4) *testl A l; 
tcirclelz = tcirclelz + coef r2 ( 1 ) *test2 A 4 + coef r2 (2) *test2 A 3 + coef r2 (3) *test2 A 2 + coefr2 
(4) *test2 A l; 

tcirclelz = tcirclelz 4- coef r3 (1) *test3 A 4 + coef r3 (2) *test3 A 3 + coef r3 (3) *test3 A 2 + coefr3 
(4) *test3 A l; 

tcirclelz = tcirclelz + coefr4 (1) *test4 A 4 + coef r4 (2 ) *test4 A 3 + coef r4 ( 3) *test4 A 2 + coefr4 
(4) *test4 A l; 

tcirclelz « tcirclelz + coef r5 (1) *test5 A 4 + coef r5 (2) *test5 A 3 + coef r5 (3) *test5 A 2 + coefrS 

(4)*test5 A l; 

testl = 2*rl; 

test2 = 0.0; 

test3 « 0.0; 

test4 = 0.0; 

test5 = 0.0; 

fprintf (fid, 1 \n\nTest of circlelz: 1 ); 

fprintf (fid, 1 \npxl = %20.12f 1 , testl) ; 

fprintf (fid, ? \npyl « %20 . 12f 1 , test2 ) ; 

fprintf (fid, f \npzl = %20. 12f 1 , test3) ; 

fprintf (fid, ' \nthl = %20 . 12f 1 , test4 ) ; 

fprintf (fid, f \nth2 = %20 . 12f 1 , tests ) ; 

fprintf (fid, 1 \ncirclelz = %20. 12f ,eval (tcirclelz) ) ; 

fprintf (fid, f \n ? ) ; 



^^yms rl thl pxl pyl pzl; 

%; dBx = diff(y2)*rlz - rly*diff (z2) /rll A 2 
%; dBy = diff(x2)*rlz - rlx*dif f (z2) /rll A 2 
%; dBz = diff (x2)*rly - rlx*dif f (y2) /rll A 2 



disp( 'finished symbolic computation of dB...'); 



%;fprintf (fid, 1 \ndBx = f ); 
%; fprintf (fid, char (dBx) ) ; 

fprintf (fid, 1 \ndBy = ' ); 
f; fprintf (fid, char (dBy) ) ; 
%; fprintf (fid, ' \ndBz = '); 
%; fprintf (fid, char (dBz) ) ; 
%; fprintf (fid, ' \n» ) ; 

%;disp( 1 finished symbolic computation of dBx, dby, dBz... 1 ); 
%; 



This next section takes too long, so it has to be converted into polynomial matrices 



pxlt = pxl; 
pylt = pyl; 
pzlt = pzl; 
thlt = thl; 
th2t = th'2; 



compute the torus equation 

- Rt = radius of torus 

- Ot = origin of torus 

- Or = origin of ring 

- Tht = angle of ring in x-y plane 

- Thtn = angle of torus in x-y plane 

- ...coordinates in torus 

- 28 points generated by for loop to 28 points generated by loop 

- looking for a 28x28 matrix to multiply each point by. . . 

- Therefor you need 28 matrices to have a solution for each row 

- Now use this matrix to generate the cumulative torus equation 

- ...with a for loop on the individual equations 



%; 
%; 
%; 

syms thcum rpt xpt ypt thpt Rt Otx Oty Otz Or Tht nr thp toruslx torusly toruslz; 
syms coef6 coef7 coef8 coef9 coeflO coefll coefl2 coefl3 coefl4 toruseqlx toruseqly 
toruseqlz; 

syms kl k2 k3 k4 tmpl Thtn thT; 
syms Rtt Otxt Otyt Otzt Thtt Thtnt; 



rl - 1.0; 
Otx = 0; 
Oty = 0; 
tz = 0; 
t = Otx + Rt*cos(Tht) ; 
ypt = Oty + Rt*sin(Tht); 
rpt = sqrt((xpt - Otx) A 2 + (ypt - Oty) A 2); 
xt = Otx + Rt*cos(Tht); 
yt = Oty + Rt*sin(Tht) ; 



%;thpt = itanl (xt, yt) ; 
syms thpt; 

thcum = thpt + Thtn; 
xO = rpt ^CTfe (thcum) ; 
^ 0 = rpt* sin (thcum) ; 

hi = Tht; 
th2 =0; y : 
nr = 8; 

thp = 2*3.14159/nr; 
Tht = thp*il; 
Rt = 20; 
toruslx = 0; 
for il = l:nr 

thpt= itanl (eval (xt) , eval (yt) ) ; 

toruslx = toruslx + circlelx; 
end 

al = [ Rt A 4 Rt A 3 Rt A 2 Rt A l;]; 
bl = [eval (toruslx) ;] ; 
Rt = 60; 
toruslx = 0; 
for il = l:nr 

thpt= itanl (eval (xt) , eval (yt) ) ; 

toruslx = toruslx + circlelx; 
end 

al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A l ; ] ; 
bl = [bl; eval (toruslx) ;] ; 
Rt = 240; 
toruslx = 0; 
for il = l:nr - 

thpt= itanl (eval (xt) , eval (yt) ) ; 

toruslx ^-faruslx + circlelx; 
end 

1 = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 2 ; ] ; 
1 - [bl; eval (toruslx) ;] ; ^ 
Rt = 2400; 
toruslx = 0; 
for il = lrnr 

thpt= itanl (eval (xt) , eval (yt) ) ; 
toruslx = toruslx + circlelx; 
end ** 
al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 2;]; 
bl = [bl; eval (toruslx) ;] ; 
%; 

The following is the first set of coefficients of the torus equation for Rt = variable 

coef6 = pinv(al)*bl; 

disp ( ' coef 6. . - 1 ) ; 

xO = Otx + Rt*cos(Tht); 

yO = Oty + Rt*sin(Tht); 

thl = Tht; 

nr = 8; 

thp = 2*3.14159/nr; 
Tht = thp*il; 
Rt = 20; 
Otx = 8*rl; 
toruslx = 0; 
for il = lrnr 

toruslx = toruslx + circlelx; 
end 

al = [Otx A 4 Otx A 3 Otx A 2 Otx A l;]; 

1 = [ (eval (toruslx) - coef 6 ( 1 ) *Rt A 4 - coef 6 (2 ) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l) ; ] ; 
^'tx = 64*rl; 
toruslx = 0; 
for il = lrnr 

toruslx = toruslx + circlelx; 
end 



al = [al; Otx A 4 Otx A 3 Otx A 2 Otx A l;]; 

bl = [bl; (eval(toruslx) - coef 6 (1) *Rt A 4 

l);]; 

Otx = -8*rl; 

•pruslx = 0; 
or il = l:nr 
toruslx = toruslx + circlelx; 
end 

al = [al; Otx A 4 Otx A 3 Otx A 2 Otx A l;]; 

bl = [bl; (eval (toruslx) - coef 6 (1) *Rt A 4 

l);]; 

Otx = -64*rl; 
toruslx = 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [al; Otx A 4 Otx A 3 Otx A 2 Otx A l;]; 

bl = [bl; (eval (toruslx) - coef 6 (1) *Rt A 4 

l);]; 



- coef6(2)*Rt A 3 - coef 6 (3) *Rt A 2 - coef6(4)*Rt A 



- coef6(2)*Rt A 3 - coef 6 (3) *Rt A 2 - coef6(4)*Rt A 



- coef6(2)*Rt A 3 - coef 6 ( 3) *Rt A 2 - coef6(4)*Rt A 



%; 

%; The following is the first set of coefficients of the torus equation for Rt = 20, Otx 

variable 

%; 

coef7 = pinv(al)*bl; 

disp ( 1 coef 7 . . . 1 ) ; 

xO = Otx + Rt*cos(Tht); 

yO = Oty + Rt*Sin(Tht); 

thl - Tht; 

nr = 8; 

thp = 2*3.14159/nr; 
Tht - thp*il; 

•t = 20; 
tx = 0; 
Oty = 8*rl; 
toruslx = 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [Oty A 4 Oty A 3 Oty A 2 Oty A l;]; 

tmp = eval (toruslx) - coef 6 (1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4) *Rt A l; 
bl = [tmp - coef7 (l)*Otx A 4 - coef 7 (2 ) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 (4 ) *Otx A l; ] ; 
Oty = 64*rl; 
toruslx = 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [al; Oty A 4 Oty A 3 Oty A 2 Oty A l;]; 

tmp = eval (toruslx) - coef 6 (1 ) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 ( 3) *Rt A 2 - coef 6 (4 ) *Rt A l; 
bl = [bl; tmp - coef 7 (1) *Otx A 4 - coef 7 (2) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 (4 ) *Otx A l; ] ; 
Oty = -8*rl; 
toruslx = 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [al; Oty A 4 Oty A 3 Oty A 2 Oty A l;]; 

tmp = eval (toruslx) - coef 6 (1) *Rt A 4 - coef 6 (2 ) *Rt A 3 - coef 6 (3 ) *Rt A 2 - coef 6 (4 ) *Rt A l; 

bl = [bl; tmp - coef 7 ( 1) *Otx A 4 - coef 7 (2) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 (4 ) *Otx A l; ] ; 

Oty = -64*rl; 

toruslx = 0; 
^^or il = l:nr 

toruslx = toruslx + circlelx; 
^end 

al = [al; Oty A 4 Oty A 3 Oty A 2 Oty A l;]; 

tmp = eval (toruslx) - coef 6 (1) *Rt A 4 - coef 6 (2 ) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6(4) *Rt A l; 
bl = [bl; tmp - coef 7 (1) *0tx A 4 - coef 7 (2) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 ( 4 ) *Otx A l; ] ; 



%; The following is the first set of coefficients of the torus equation for Rt = 20, Otx = 
0, Oty = variable 

ioef8 = pinv(al)*bl; 

disp( 1 coef8. . . f ) ; 

xO = Otx + Rt*cos(Tht); 

yO = Oty + Rt*sin(Tht); 

thl = Tht; 

nr = 8; 

thp = 2*3.14159/nr; 
Tht = thp*il; 
Rt = 20; 
Otx = 0; 
Oty = 0; 
otz = 8*rl; 
toruslx = 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [Otz A 4 Otz A 3 Otz A 2 Otz A l;]; 

tmp » eval (toruslx) - coef 6 ( 1) *Rt A 4 - coef 6 (2 ) *Rt A 3 - coef 6 ( 3) *Rt A 2 - coef 6 (4 ) *Rt A l; 
tmp = tmp~- coef7 (1) *Otx A 4 - coef 7 (2) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 (4 ) *Otx A l; 
bl = [tmp - coef8(l)*Oty A 4 - coef 8 (2) *Oty A 3 - coef 8 (3) *Oty A 2 - coef 8 (4 ) *Oty A l; ] ; 
Otz = 64*rl; 
toruslx = 0; 
for il = l:nr 

torus 1x7^ toruslx + circlelx; 
end 

al — [al; Otz A 4 Otz A 3 Otz A 2 Otz A l;]; 

tmp = eval(toruslx) - coef 6 (1) *Rt A 4 ' - coef 6 (2 ) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l; 
^mp = tmp - coef7 (1) *Otx A 4 - coef 7 (2) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 (4 ) *Otx A l ; 
Pi = [bl; tmp - coef8(l)*Oty A 4 - coef 8 (2 ) *Oty A 3 - coef 8 (3) *Oty A 2 - coef 8 ( 4 ) *Oty A l; ] ; 
Otz = -8*rl; 
toruslx = 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [al; Otz A 4 Otz A 3 Otz A 2 Otz A l;]; 

tmp = eval(toruslx) - coef 6 (1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l ; 

tmp = tmp - coef7 (1) *Otx A 4 - coef 7 (2) *Otx A 3 - coef 7 ( 3) *Otx A 2 - coef 7 (4 ) *Otx A l; 

bl = [bl; tmp - coef8(l)*Oty A 4 - coef 8(2) *Oty A 3 - coef 8 (3) *Oty A 2 - coef 8 ( 4 ) *Oty A l; ] ; 

Otz = -64*rl; 

toruslx = 0; 

for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [al; Otz A 4 Otz A 3 Otz A 2 Otz A l;] ; 

tmp = eval(toruslx) - coef 6 (1) *RtV - coef 6 (2 ) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l; 

tmp = tmp - coef7 (1) *Otx A 4 - coef 7 (2) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 ( 4 ) *Otx A l ; 

bl = [bl; tmp - coef 8 ( 1 ) *Oty A 4 - coef 8 (2 ) *Oty A 3 - coef 8 (3) *Oty A 2 - coef8 (4 ) *Oty A l; ] ; 



%; 

%; The following is the first set of coefficients of the torus equation for Rt = 20 f Otx = 

0, Oty = 0, Otz = var. 

%; 

coef 9 = pinv(al)*bl; 
disp( 1 coef 9. . . 1 ) ; 
%; 

x0 = Otx + Rt*cos(Tht); 

0 = Oty + Rt*sin(Tht); 
thl = Tht; 
nr = 8; 

thp - 2*3.14159/nr; 
Tht = thp*il; 



kl - 10; _ „: 
Jc2 = 2; 
k3 = 2; 
k4 = -2; 

•t = kl*rl; | 
tx = k2*rl; * 
Oty = k3*rl; 
Otz = k4*rl; 
pxl = 2*rl; 
toruslx = 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

tmpl = eval (toruslx) - coef 6 (1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l; 
tmpl » tmpl - coef7 (1) *Otx A 4 - coef 7 (2 ) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 (4 ) *Otx A l; 
tmpl = tmpl - coef 8 (1) *Oty A 4 - coef 8 (2) *Oty A 3 - coef 8 (3) *Oty A 2 - coef 8 (4 ) *Oty A l; 
tmpl = tmpl - coef9(l) *Otz A 4 - coef 9 (2) *Otz A 3 - coef 9 (3) *Otz A 2 - coef 9 ( 4 ) *Ot z A l ; 



bl = [tmpl; ] ; 
pxl = 8*rl; 
toruslx = 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [al; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 
bl = [bl; tmpl;] ; 
pxl = 32*rl; 
toruslx ^*Vr 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

1 = [al; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

1 = [bl; tmpl;]; 
pxl = 200*rl; 
toruslx - 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [al; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

bl = [bl; tmpl;] ; 

%; 

%; coefficients of the torus equation for (...) pxl = variable 
%; 

coeflO = pinv(al)*bl; 
disp( 'coef 10. . . 1 ) ; 



pxl = 2*rl; 
pyl = 2*rl; 
toruslx = 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

bl = [tmpl - coeflO(l)*pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; ] ; 
pyl = 8*rl; 
toruslx = 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

= [al; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 
il = [bl; tmpl - coef 10(1) *pxl A 4 - coef 10 (2 ) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; ] ; 
pyl = 32*rl; 
toruslx = 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 



end 

al = [al; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

bl = [bl; tmpl - coef 10 (1) *pxl A 4 - coef 10 (2 ) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; ] ; 
pyl = 200*rl; 
toruslx = 0; 
tor il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [al; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

bl = [bl; tmpl - coef 10 (1) *pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; ] ; 
%; 

%; coefficients of the torus equation for (...) pyl = variable 
%; 

coefll = pinv(al)*bl; 
disp( 'coefll. . . ' ) ; 
pxl = 2*rl; 
pyl = 2*rl; 
pzl = 2*rl; 
toruslx = 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

tmp = tmpl - coeflO (l)*pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 
bl = [tmp - coefll (1) *pyl A 4 - coef 11 (2 ) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 ( 4 ) *pyl A l; ] ; 
pzl = 8*rl; 
toruslx = 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [al; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

tmp = tmpl - coef 10(1) *pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 
)1 = [bl; tmp - coefll (l)*pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4 ) *pyl A l; ] ; 
Izl = 32*rl; 
'toruslx = 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [al; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

tmp = tmpl - coefl0(l)*pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 

bl = [bl; tmp - coef 11 (1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 ( 3) *pyl A 2 - coef 11 (4 ) *pyl A l; ] ; 

pzl - 200*rl; 

toruslx = 0; 

for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [al; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

tmp = tmpl - coef 10(1) *pxl A 4 - coef 10 ( 2) *pxl A 3 - coef 10 (3) *pxl A 2 - coeflO (4) *pxl A l; 

bl = [bl; tmp - coef 11 (1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4) *pyl A l; ] ; 

coefficients of the torus equation for (...) pzl = variable 

%; 

coefl2 = pinv(al)*bl; 
disp( 'coef 12. . . 1 ) ; 
pxl = 2*rl; 
pyl = 2*rl; 
pzl = 2*rl; 
Tht = 3.14159/2; 
toruslx = 0; 
for il = l:nr 
toruslx = toruslx + circlelx; 
"end 

al = [Tht A 4 Tht A 3 Tht A 2 Tht A l;]; 

tmp = tmpl - coefl0(l)*pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 ( 3) *pxl A 2 - coeflO (4) *pxl A l; 
tmp = tmp - coefll (l)*pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 ( 4 ) *pyl A l ; 



bl = [tmp - coefl2(l)*pzl A 4- coef 12 (2) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 (4 ) *pzl A l; ] ; 
Tht = 3.14159; 
toruslx = 0; 
for il ^WOxr" 

toruslx = toruslx + circlelx; 
Ind 

al = [al; > Tht A 4 Tht^3 Tht A 2 Tht A l;]; 

tmp = tmpi - coeflO(l)*pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 

tmp = tmp - coef 11 (1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4) *pyl A l; 

bl = [bl; tmp - coef 12 (1) *pzl A 4- coefl2 (2) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 (4 ) *pzl A l; ] ; 

Tht = 3*3.14159/2; 

toruslx = 0; 

for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [al; Tht A 4 Tht A 3 Tht A 2 Tht A l;]; 

tmp = tmpl - coefl0(l)*pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4) *pxl A l; 

tmp = tmp - coefll (1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4) *pyl A l; 

bl « [bl; tmp - coef 12 (1) *pzl A 4- coef 12 (2) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 (4 ) *pzl A l ; ] ; 

Tht = 2*3.14159; 

toruslx = 0; 

for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [al; Tht A 4 Tht A 3 Tht A 2 Tht A l;]; 

tmp = tmpl - coeflO (1) *pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 ( 3) *pxl A 2 - coef 10 (4 ) *pxl A l; 

tmp = tmp - coefll (l)*pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4 ) *pyl A l; 

bl = [bl; tmp - coef 12 (1) *pzl A 4- coef 12 (2) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 (4 ) *pzl A l; ] ; 



coefficients of the torusr equation for (...) Tht = variable 



^oefl3 = pinv(al)*bl; 

Usp( 1 coef 13. . . 1 ) ; 
"pxl = 2*rl; 
pyl = 2*rl; 
pzl = 2*rl; 
Tht = 0.0; 
Thtn = 3.14159/2; 
toruslx = 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [Thtn A 4 Thtn A 3 Thtn A 2 Thtn A l;]; 

tmp = tmpl - coeflO(l)*pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 ( 3) *pxl A 2 - coef 10 (4 ) *pxl A l; 
tmp = tmp - coefll (l)*pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4 ) *pyl A l; 
tmp = tmp - coefl2(l)*pzl A 4- coef 12 (2) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 (4 ) *pzl A l; 
bl = [tmp - coefl3(l)*Tht A 4 - coef 13 (2) *Tht A 3 - coef 13 (3) *Tht A 2 - coef 13 (4 ) *Tht A l; ] ; 
Thtn = 3.14159; 
toruslx = 0; 
for il = l:nr 

toruslx = toruslx + circlelx; 
end 

al = [al; Thtn A 4 Thtn A 3 Thtn A 2 Thtn A l;]; 

tmp = tmpl - coefl0(l)*pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 
tmp = tmp - coefll (1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4 ) *pyl A l; 
tmp = tmp - coef 12 (1) *pzl A 4- coef 12 (2 ) *pzl A 3 - coef 12 ( 3) *pzl A 2 - coef 12 (4) *pzl A l; 
bl = [bl; tmp - coef 13 (1) *Tht A 4 - coef 13 (2) *Tht A 3 - coef 13 (3) *Tht A 2 - coef 13 (4 ) *Tht A l; ] 
Thtn - 3*3.14159/2; 
toruslx = 0; 
for il = l:nr 
toruslx = toruslx + circlelx; 
"end 

al = [al; Thtn A 4 Thtn A 3 Thtn A 2 Thtn A l;]; 

tmp = tmpl - coeflO(l)*pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 ( 3) *pxl A 2 - coeflO (4) *pxl A l; 
tmp = tmp ~ coefll (l)*pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 ( 3) *pyl A 2 - coef 11 (4 ) *pyl A l; 



tmp = tmp - coefl2(l)*pzl A 4- coef 12 (2) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 (4 ) *pzl A l; 
bl = [bl; tmp - coef 13 (1) *Tht A 4 - coef 13 (2 ) *Tht A 3 - coef 13 ( 3) *Tht A 2 - coef 13 (4 ) *Tht A l; ] ; 
Thtn - 2*3.14159; 
toruslx = 0; 
or il = l:nr 
toruslx = toruslx + circlelx; 
end 

al = [al; Thtn A 4 Thtn A 3 Thtn A 2 Thtn A l;]; 

tmp = tmpl - coef 10(1) *pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coeflO (4) *pxl A l; 
tmp = tmp - coefll (1) *pyl A 4 - coef 11 (2 ) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 ( 4 ) *pyl A l; 
tmp = tmp - coefl2(l)*pzl A 4- coef 12 (2) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 ( 4 ) *pzl A l; 
bl = [bl; tmp - coef 13 (1) *Tht A 4 - coef 13 (2) *Tht A 3 - coef 13 (3) *Tht A 2 - coef 13 (4) *Tht A l; ] ; 



coefficients of the torus equation for (...) Thtn = variable 



% 
% 
% 

coefl4 = pinv(al)*bl; 
disp( 1 coef 14. . . ' ) ; 



Now, finally the complete torus equation 



coef 6 = eval (coef 6); 
coef7 = eval(coef7); 
coef 8 = eval (coef 8); 
coef 9 = eval (coef 9); 
coeflO = eval(coeflO) 
coefll = eval (coefll) 
coef 12 = eval (coef 12) 
coef 13 = eval (coef 13) 
coef 14 = eval (coef 14) 

yms Rtt Otxt Otyt Otzt pxlt pylt pzlt Thtt Thtnt; 

toruseqlx = coef 6 (1) *Rtt A 4 + coef 6 (2 ) *Rtt A 3 + coef 6 ( 3) *Rtt A 2 + coef 6 (4 ) *Rtt A l; 
toruseqlx = toruseqlx + coef 7 (1) *Otxt A 4 + coef7 (2) *Otxt A 3 + coef 7 (3) *Otxt A 2 + coef7(4) 
*Otxt A l; 

toruseqlx = toruseqlx + coef8 (1) *Otyt A 4 + coef 8 (2) *Otyt A 3 + coef 8 (3) *Otyt A 2 + coef8(4) 
*Otyt A l; 

toruseqlx = toruseqlx + coef 9 ( 1) *Otzt A 4 + coef 9 (2) *Otzt A 3 + coef 9 (3) *Otzt A 2 + coef 9 (4) 
*Otzt A l; 

toruseqlx = toruseqlx + coef 10 (1 ) *pxlt A 4 + coef 10 (2 ) *pxlt A 3 + coef 10 (3) *pxlt A 2 + coeflO (4) 
*pxlt A l; 

toruseqlx = toruseqlx + coef 11 ( 1 ) *pylt A 4 + coef 11 (2) *pylt A 3 + coef 11 (3) *pylt A 2 + coefll (4) 
*pylt A l; 

toruseqlx = toruseqlx + coef 12 ( 1) *pzlt A 4 + coef 12 (2) *pzlt A 3 + coef 12 (3) *pzlt A 2 + coefl2(4) 
*pzlt A l; 

toruseqlx = toruseqlx + coef 13 (1) *Thtt A 4 + coef 13 (2) *Thtt A 3 + coef 13 (3) *Thtt A 2 + coef 13 (4) 
*Thtt A l; 

toruseqlx = toruseqlx + coef 14 ( 1 ) *Thtnt A 4 + coef 14 (2) *Thtnt A 3 + coef 14 (3) *Thtnt A 2 + coefl4 
(4)*Thtnt A l; 

disp ( 1 finished computation of toruseqlx... 1 ); 
fprintf (fid, ' \n\ntoruseqlx = Sum of: \n f ); 
for il = 1:4 

fprintf (fid, ? %20 . 12f 1 , coef 6 (il) ) ; 

fprintf (fid, f Rtt A %li 1 , il ) ; 

fprintf (fid, 1 \n f ) ; 
end 

for il = 1:4 

fprintf (fid, 1 %20 . 12f 1 , coef 7 (il) ) ; 

fprintf (fid, f Otxt A %li 1 , il) ; 

fprintf (fid, '\n' ) ; 
end 

for il = 1:4 . _ 

fprintf (fid, I %20.12f i : ;> Qoef 8 (il) ) ; .^-.^ 



fprintf (fid, ? Otyt A %li 1 , il) ; 
f print f (fid, *\n' ) ; 
end 

for il = 1:4 

fprintf (fid, 1 %20. 12f 1 , coef 9 (il) ) ; 

fprintf (fid, 1 Otzt A %li ' , il) ; 

fprintf (fid, 'Xn' ) ; 
end 

for il = 1:4 

fprintf (fid, 1 %20 . 12f 1 , coef 10 (il) ) ; 

fprintf (fid, 1 pxlt A %li ' , il ) ; 

fprintf (fid, 'Xn') ; 
end 

for il = 1:4 

fprintf (fid, '%20.12f 1 , coef 11 (il) ) ; 

fprintf (fid, 1 pylt A %li • , il) ; 

fprintf (fid, 1 \n f ) ; 
end 

for il = 1:4 

fprintf (fid, ' %20 . 12f 1 , coef 12 (il) ) ; 

fprintf (fid, 1 pzlt A %li ' , il ) ; 

fprintf (fid, '\n'); 
end 

for il = 1:4 

fprintf (fid, 1 %20 . 12f f , coef 13 (il ) ) ; 

fprintf (fid, 1 Thtt A %li 1 , il) ; 

fprintf (fid, 'Nn') ; 
end 

for il = 1:4 

fprintf (fid, 1 %20 . 12f 1 , coef 14 (il) ) ; 

fprintf (fid, 1 Thtnt A %li 1 , -il) ; 

fprintf (fid, 1 \n») ; 
|nd 

fprintf (fid, ' \n\nTest of toruslx'); 

Rtt = 5.0; 

Otxt =0.0; 

Otyt = O.O7 

Otzt = 0.0; 

pxlt = 0.0; 

pyit- = a*rtr? ;r 

pzlt = 0.0; 
Thtt = 0.0; 
Thtnt = 0.0; 

fprintf (fid, ? \nRtt = %20 . 12f f , Rtt ) ; 

fprintf (fid, ' \nOtx = %20 . 12f 1 , Otxt) ; 

fprintf (fid, '\nOty = %20. 12f 1 ,Otyt) ; 

fprintf (fid, '\nOtz = %20 . 12f f , Otzt ) ; 

fprintf (fid, ! \npxlt = %20.12f 1 ,pxlt) ; 

fprintf (fid, 1 \npylt = %20 . 12f f , pylt ) ; 

fprintf (fid, f \npzlt = %20 . 12f ' , pzlt ) ; 

fprintf (fid, 1 \nThtt = %20 . 12f 1 , Thtt ) ; 

fprintf (fid, 1 \nThtnt = %20 . 12f 1 , Thtnt ) ; 

fprintf (fid, f \ntoruslx = %20 . 12f ? , eval (toruslx) ) ; 

fprintf (fid, 'Xn 1 ) ; 

fprintf (fid, 1 \n\nTest of toruslx 1 ); 
Rtt = 5.0; 
Otxt = 0.0; 
Otyt =0.0; 
Otzt = 0.0; 
xlt = 10.0; 
ylt = 0.0; 
pzlt = 0.0; 
Thtt =0.0; 
Thtnt = 0.0; 

fprintf (fid, 1 XnRtt = %20 . 12f 1 , Rtt ) ; 



m 



fprintf (fid, '\nOtx = %20 . 12f ' , Otxt ) ; 
fprintf (fid, 1 \nOty = %20 . 12f ' , Otyt ) ; 
fprintf (fid, f ,\n6tz = %20 . 12f 1 , Otzt ) ; 
fprintf (fid, »\npxlt = %20. 12f ' ,pxlt) ; 

printf (fid, f \npylt = %20.12f ' ,pylt) ; 

printf (fid, f \npzlt = %20. 12f ' ,pzlt) ; 
fprintf (fid, 1 \nThtt = %20. 12f f , Thtt) ; 
fprintf (fid, '\nThtnt = %20. 12f\ Thtnt) ; 
fprintf (fid, 1 \ntoruslx = %20 . 12f 1 , eval ( toruslx) ) ; 
fprintf (fid, f \n f ) ; 

fprintf (fid, 1 \n\nTest of toruslx 1 ); 

Rtt = 5.0; 

Otxt = 0.0; 

Otyt = 0.0; 

Otzt = 0.0; 

pxlt = 0.0; 

pylt = 0.0; 

pzlt = 10.0; 

Thtt = 0.0; 

Thtnt = 0.0; 

fprintf (fid, ' \nRtt = %20 . 12f ' , Rtt ) ; 

fprintf (fid, 1 \nOtx = %20 . 12f ' , Otxt ) ; 

fprintf (fid, 1 \nOty - %20 . 12f 1 , Otyt ) ; 

fprintf (fid, f \nOtz = %20 . 12f ' , Otzt ) ; 

fprintf (fid, '\npxlt = %20 . I2f 1 , pxlt ) ; 

fprintf (fid, f \npylt = %20 . 12f 1 , pylt ) ; 

fprintf (fid, ? \npzlt = %20. 12f 1 , pzlt) ; 

fprintrftfTO) 'AnThtt = %20. TTf ' , Thtt ) ; . 

fprintf (fid, f \nThtnt = %20. 12f f , Thtnt) ; 

fprintf (fid, 1 \ntoruslx = %20. 12f , eval (toruslx) ) ; 

fprintf (fid, '\n % ) ; 

fprintf (fid, 1 \n\nTest of toruslx 1 ); 
.tt = 5.0; 
txt = 0.0; 
Otyt = 0.0; 
Otzt = 0.0; 
pxlt = -10.0; 
pylt = 0.0; 
pzlt = 0.0; 
Thtt = 0.0;. 
Thtnt « 0.0; 

fprintf (fid, ? \nRtt = %20 . 12f 1 , Rtt ) ; 

fprintf (fid, 1 \nOtx = %20 . 12f 1 , Otxt ) ; 

fprintf (fid, 1 \nOty = %20.12f \Otyt) ; 

fprintf (fid, '\nOtz = %20 . 12f 1 , Otzt ) ; 

fprintf (fid, 1 \npxlt = %20. 12f Vpxlt ) ; 

fprintf (fid, 1 \npylt = %20 . 12f 1 , pylt ) ; 

fprintf (fid, '\npzlt = %20. 12f 1 ,pzlt) ; 

fprintf (fid, f \nThtt = %20. 12f » , Thtt ) ; 

fprintf (fid, ' \nThtnt = %20 . 12f Thtnt ) ; 

fprintf (fid, f \ntoruslx = %20 . 12f , eval (toruslx) ) ; 

fprintf (fid,. ! \n f ) ;. 

fprintf (fid, 1 \n\nTest of toruslx 1 ); 

Rtt = 5.0; 

Otxt = 0.0; 

Otyt = 0.0; 

Otzt = 0.0; 

pxlt = 0.0; 

pylt = 0.0; 

pzlt = -10.0; 

Thtt = 0.0; 

htnt = 0.0; 
fprintf (fid, '\nRtt = %20 . 12f T , Rtt ) ; 
fprintf (fid, ? \nOtx = %20-.12f ' , Otxt) ;- 
fprintf (fid, '\nOty * %20 . 12f 1 , Otyt ) ; 
fprintf (fid, f \nOtz = %20 . 12f 1 , Otzt ) ; 



fprintf (fid, f \npxlt = %20 . 12f 1 , pxlt ) ; 
fprintf (fid, 1 \npylt = %20 . 12f f , pylt ) ; 
fprintf. (fid, 1 \npzlt = %20.12f 1 , pzlt) ; 
fprintf (fid, 'Vfthtt = %20 . 12f 1 , Thtt ) ; 

•printf (fid, 1 \nThtnt = %20. 12f 1 , Thtnt) ; 
printf (fid, f \ntoruslx = %20 . 12f 1 , eval ( toruslx) ) ; 
fprintf (fid,. 1 \rr' ) ;. 



syms coef6 coef7 coef8 coef9 coeflO coefll coef!2 coef!3 coefl4; 



rl = 1.0;. 
Otx = 0; 
Oty = 0; 
Otz = 0; 

xpt = Otx + Rt*cos (Tht) ; 
ypt = Oty + Rt*sin(Tht); 

rpt = sqrt((xpt - Otx) A 2 + (ypt - Oty) A 2); 
xt = xpt; 
yt = ypt; 

%; thpt = itanl(xt,yt) ; 
syms thpt; 

thcum = thpt + Thtn; 
xO = rpt*cos (thcum) ; 
yO = rpt*sin (thcum) ; 
thl = Tht; 
th2 = 0; 
nr = 8; 

thp = 2*3.14159/nr; 
Tht = thp.*il; 
torusly = 0; 

• t " 20; 
or il = l:nr 

thpt = itanl (xt, yt) ; 

torusly = torusly + circlely; 

end 

al = [ Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 

bl = [eval (toruslx) ;] ; 

Rt - 60; 

for il = l:nr 

thpt = itanl (xt, yt) ; 

torusly = torusly + circlely; 
end 

al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 
bl = [bl; eval (torusly) ;] ; 
Rt = 240; 
for il = l:nr 

thpt = itanl (xt, yt) ; 

torusly = torusly + circlely; 
end 

al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 2 ; ] ; 
bl = [bl; eval (torusly) ;] ; 
Rt = 2400; 
for il = l:nr 

thpt = itanl (xt, yt) ; 

torusly = torusly + circlely; 
end 

al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 2 ; ] ; 
bl = [bl; eval (torusly) ;] ; 

•V 
|fe; The following is the first set of coefficients of the torus equation for Rt = variable 
% ; 

coef6 = pinv(al)*bl; 

disp ( 1 coef 6 . . . 1 ) ; 

xO = Otx + Rt*cos(Tht) — 



yO = Oty + Rt*sin(Tht); 
thl = Tht; 
nr = 8; 

thp = 2*3.14159/nr; 



• V; 



Otx = 8*rl; 




torusly = 0; 
for il = l:nr 

torusly = torusly + circlely; 
end 

al = [Otx A 4 Otx A 3 Otx A 2 Otx A l;]; 

bl - [(eval (torusly) - coef 6 (1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l) ; ] ; 
Otx = 64*rl; 
torusly = 0; 
for il = l:nr 

torusly = torusly + circlely; 
end 

al = [al; Otx A 4 Otx A 3 Otx A 2 Otx A l;]; 

bl = [bl; (eval (torusly) - coef 6 ( 1 ) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef6(4)*Rt A 
l);]; 

Otx = -8*rl; 
torusly = 0; 
for il = l:nr 

torusly^= torusly + circlely; 
end 

al = [al; Otx A 4 Otx A 3 Otx A 2 Otx A l;]; 

bl = [bl; (eval (torusly) - coef 6 (1 ) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef6(4)*Rt A 
l);]; 

Otx = -64*rl; 
torusly = 0; 
for il = l:nr 

torusly = torusly + circlely; 
nd 

al = [al; Otx A 4 Otx A 3 Otx A 2 Otx A l;]; 

bl = [bl; (eval (torusly) - coef 6 (1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef6(4)*Rt A 
D;3; 

%; 

%; The following is the first set of coefficients of the torus equation for Rt = 20, Otx = 

variable 

%; 

coef7 = pinv(al)*bl; 

disp ( f coef 7 . . . ? ) ; 

xO - Otx + Rt*cos(Tht); 

yO = Oty + Rt*sin(Tht) ; 

thl = Tht; 

nr = 8; 

thp = 2*3.14159/nr; 

Tht = thp*il; 

Rt = 20; 

Otx = 0; 

Oty = 8*rl; 

torusly = 0; 

for il = l:nr 

torusly = torusly + circlely; 
end 

al = [Oty A 4 Oty A 3 Oty A 2 Oty A l;]; 

txnp = eval (torusly) - coef 6 ( 1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l; 
bl = [tmp - coef7 (l)*Otx A 4 - coef 7 (2 ) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 (4 ) *Otx A l; ] ; 




torusly- = torusly + circlely; 
end 



al = [al; Oty A 4 Oty A 3 Oty A 2 Oty A l;]; 



tmp = eval (torusly) - coef 6 .( 1 ) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l; 
bl = [bl; tmp - coef7 (1) *0tx A 4 - coef 7 (2 ) *0tx A 3 - coef 7 (3) *0tx A 2 - coef7 (4) *Otx A l; ] ; 
Oty = -8*rl; 
torus ly.j*?&r 
or il = l:nr 
torusly = torusly + circlely; 
end ^. ; ~ 

al = [al;'~0ty A 4 0ty A 3 0ty A 2 Oty A l;]; 

tmp = eval (torusly) - coef 6 (1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l; 
bl = [bl; tmp - coef 7 (1) *0tx A 4 - coef 7 (2) *0tx A 3 - coef 7 (3) *0tx A 2 - coef 7 (4 ) *Otx A l; ] ; 
Oty = -64*rl; 
torusly = 0; 
for il = l:nr 

torusly = torusly + circlely; 
end 

al = [al; Oty A 4 Oty A 3 Oty A 2 Oty A l;]; 

tmp = eval (torusly) - coef 6 (1) *Rt A 4 - coef 6 (2 ) *Rt A 3 - coef 6 ( 3) *Rt A 2 - coef 6 ( 4 ) *Rt A l; 
bl = [bl; tmp - coef 7 ( 1) *Otx A 4 - coef 7 (2 ) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 ( 4 ) *Otx A l; ] ; 

%; 

%; The following is the first set of coefficients of the torus equation for Rt = 20, Otx 

0, Oty = variable 

%; 

coef 8 = pinv(al)*bl; 

disp ( 1 coef 8 ...'); 

xO = Otx + Rt*cos (Tht) ; 

yO = Oty + Rt*sin(Tht); 

thl - Tht; 

nr = 8; 

thp = 2*3.14159/nr; 

Tht = thp*il; 

Rt = 20; 
^^)tx = 0; 
^ty = 0; 

Otz = 8*rl; 

torusly = 0; 

for il = l:nr 

torusly = torusly + circlely; 

end 

al = [Otz A 4 Otz A 3 Otz A 2 Otz A l;]; 

tmp = eval (torusly) - coef 6 (1) *Rt A 4 - coef € (2 ) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4) *Rt A l; 
tmp = tmp - coef7 (1) *Otx A 4 - coef 7 (2 ) *Otx A 3 - coef 7 ( 3 ) *Otx A 2 - coef 7 ( 4 ) *Otx A l; 
bl = [tmp - coef8 (1) *Oty A 4 - coef 8 (2 ) *Oty A 3 - coef 8 (3) *Oty A 2 - coef 8 ( 4 ) *Oty A l; ] ; 
Otz = 64*rl; 
torusly = 0; 
for il = l:nr 

torusly = torusly + circlely; 
end 

al = [al; Otz A 4 Otz A 3 Otz A 2 Otz A l;]; 

tmp = eval(torusly) - coef 6 (1) *Rt A 4 - coef 6 (2 ) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 ( 4 ) *Rt A l ; 

tmp = tmp - coef7 (1) *Otx A 4 - coef 7 (2) *Otx A 3 - coef 7 ( 3) *Otx A 2 - coef 7 (4 ) *Otx A l; 

bl = [bl; tmp - coef 8 (1) *Oty A 4 - coef 8 (2) *Oty A 3 - coef 8 (3) *Oty A 2 - coef 8 (4) *Oty A l; ] ; 

Otz = -8*rl; 

torusly = 0; 

for il = l:nr 

torusly = torusly + circlely; 
end 

al = [al; Otz A 4 Otz A 3 Otz A 2 Otz A l;]; 

tmp = eval (torusly) - coef 6 (1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l; 
tmp = tmp - coef7(l)*Otx A 4 - coef 7 (2 ) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 (4 ) *Otx A l; 

•bl = [bl; tmp - coef 8 (1) *Oty A 4 - coef 8 (2) *Oty A 3 - coef 8 (3) *Oty A 2 - coef 8 (4 ) *Oty A l; ] ; 
ptz = -64*rl; 
torusly = 0; 
for il = l:nr 

torusly = torusly + circlely; 
end 



al = [al; Otz A 4 Otz A 3 Otz A 2 Otz A l;]; 

tmp = eval (torusly) - coef 6 (1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 ( 4 ) *Rt A l; 

tmp = tmp - coef7 (1) *0tx A 4 - coef? (2) *0tx A 3 - coef 7 (3) *Otx A 2 - coef 7 ( 4 ) *Otx A l; 

bl = [bl; tmp - coef 8 (1) *0ty A 4 - coef 8 (2) *0ty A 3 - coef 8 (3) *0ty A 2 - coef 8 ( 4 ) *Oty A l; ] ; 



%; The following is the first set of coefficients of the torus equation for Rt = 20, Otx = 
0, Oty = 0, Otz = var. 

%; 

coef9 - ginv(al) *bl; 
disp( 'coef 9. . . ' ) ; 
%; 

x0 = Otx + Rt*cos(Tht); 
yO = Oty + Rt*sin(Tht); 
thl = Tht; 
nr = 8; 

thp = 2*3.14159/nr; 

Tht = thp*il; 

Rt = kl*rl; 

Otx = k2*rl; 

Oty = k3*rl; 

Otz = k4*rl; 

pxl = 2*rl; 

torusly = 0; 

for il = l:nr 

torusly = torusly + circlely; 
end 

al = [pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

tmpl = eval (torusly) - coef 6 (1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l; 
tmpl = tmpl - coef7 (1) *Otx A 4 - coef 7 (2 ) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 ( 4 ) *Otx A l; 
tmpl = tmpl - coef8 (1) *Oty A 4 - coef 8 (2) *Oty A 3 - coef 8 (3) *Oty A 2 - coef8 (4) *Oty A l; 
tmpl = tmpl - coef 9(1) *Otz A 4 - coef 9 (2 ) *0tz A 3 - coef 9 (3) *Otz A 2 - coef 9 (4) *Otz A l; 



1 - [tmpl;] ; 
pxl = 8*rl; 
torusly = 0; 
for il = l:nr 

torusly = torusly + circlely; 
end 



al = [al; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 
bl = [bl; tmpl; ] ; 
pxl = 32*rl; 
torusly = 0; 
for il = l:nr 

torusly = torusly + circlely; 
end 

al = [al; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 
bl = [bl; tmpl; ] ; 
pxl = 200*rl; 
torusly = 0; 
for il = l:nr 

torusly = torusly + circlely; 
end 

al = [al; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

bl = [bl; tmpl;] ; 

%; 

%; coefficients of the torus equation for (...) pxl = variable 
%; 

coeflO = pinv(al)*bl; 
disp ( 1 coef 10 ...'); 

pxl = 2*rl; 
pyl = 2*rl; 
torusly = 0; 
for il = l:nr 



torusly_= torusly + circlely; 
end 

al = [pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

bl - [trrtgj - cdfrfl0(11*pxl^f - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; ] / 
h yl = 8*rl; 

orusly = 0; 
for il = l:nr 

torusly = torusly + circlely; 
end 

al = [al;pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

bl = [bl; tmpl - coeflO(l)*pxl A 4 - coef 10 (2 ) *pxl A 3 - coef 10 ( 3) *pxl A 2 - coef 10 (4 ) *pxl A l; ] ; 
pyl = 32*rl; 
torusly = 0; 
for il = l:nr 

torusly = torusly + circlely; 
end 

al - [al; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

bl = [bl; tmpl - coef 10 ( 1 ) *pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4) *pxl A l; ] ; 
pyl = 200*rl; 
torusly = 0; 
for il = l:nr 

torusly = torusly + circlely; 
end 

al = [al; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

bl = [bl; tmpl - coef 10 (1) *pxl A 4 - coef 10 (2 ) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4) *pxl A l; ] ; 
ft; 

ft; coefficients of the torus equation for (...) pyl = variable 
ft; 

coef 11 = pinv(al)*bl; 
disp ( ' coef 11 . . . 1 ) ; 
pxl = 2*rfr 
pyl = 2*rl; 

>zl = 2*rl; 

torusly = 0; 
for il = l:nr 

torusly = torusly + circlely; 
end 

al = [pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

tmp = tmpl - coef 10 (1) *pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 
bl « [tmp - coefll(l)*pyl A 4 - coef 11 (2 ) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 ( 4 ) *pyl A l; ] ; 
pzl = 8*rl; 
torusly = 0; 
for il = l:nr 

torusly = torusly + circlely; 
end 

al = [al; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

tmp = tmpl - coefl0<l)*pxl A 4 - coef 10 (2 ) *pxl A 3 - coef 10 (3) *pxl A 2 - coeflO (4) *pxl A l; 

bl = [bl; tmp - coef 11 ( 1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4 ) *pyl A l; ] ; 

pzl = 32*rl; 

torusly = 0; 

for il = l:nr 

torusly = torusly + circlely; 
end 

al - [al; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

tmp = tmpl - coeflO (1) *pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coeflO (4) *pxl A l; 

bl = [bl; tmp - coef 11 (1) *pyl A 4 - coef 11 ( 2 ) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 ( 4 ) *pyl A l; ] ; 

pzl = 200*rl; 

torusly = 0; 

for il = l:nr 

torusly = torusly + circlely; 
md 

)l = [al; pzl A 4 pzl A 3 pzl A 2 pzl A l;j; 
'fcmp = tmpl - coef 10 (1) *pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 
bl - [bl; tmp - coef 11 (1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4 ) *pyl A l; ] ; 



ft; 



%; coefficients of the torus equation for 
%; 



, ) pzl = variable 



coefl2 = pinv(al)*bl; 
disp( 'coefl2. . . 1 ) ; 
^|xl = 2*rl; 
5yl = 2*rl; 
pzl = 2*rl; 
Tht = 3.14159/2; 
torusly = 0; 
for il = l:nr 

torusly = torusly + circlely; 
end 

al = [Tht A 4 Tht A 3 Tht A 2 Tht A l;]; 

tmp = tmpl - coeflO (1) *pxl A 4 - coef 10 (2 ) *pxl A 3 

tmp = tmp - coef 11 (1) *pyl A 4 - coef 11 (2) *pyl A 3 - 

bl = [tmp.- coefl2(l)*pzl A 4- coef 12 (2 ) *pzl A 3 - 

Tht = 3.14159; 

torusly = 0; 

for il = l:nr 

torusly = torusly + circlely; 
end 

al = [al; Tht A 4 Tht A 3 Tht A 2 Tht A l;]; 

tmp = tmpl - coef 10 (1) *pxl A 4 - coef 10 (2) *pxl A 3 

tmp = tmp - coef 11(1) *pyl A 4 - coef 11 (2) *pyl A 3 - 

bl = [bl; tmp - coef 12 ( 1 ) *pzl A 4- coef 12 (2) *pzl A 

Tht = 3*3.14159/2; 

torusly = 0; 

for il = -Tftnr 

torusly = torusly + circlely; 
end 

al = [al;.^Tht A 4 y Tht A 3 Tht A Z : Tht A l; ] ; 
tmp = tmpl - coef 10 (1) *pxl A 4 - coef 10 (2) *pxl A 3 
tmp = tmp - coef 11 (1) *pyl A 4 - coef 11 (2) *pyl A 3 
|l = [bl; tmp - coefl2(l)*pzl A 4- coef 12 (2 ) *pzl A 
^Tht = 2*3.14159; 
torusly = 0; 
for il = l:nr 

torusly - torusly + circlely; 
end 

al = [al; Tht A 4 Tht A 3 Tht A 2 Tht A l;]; 
tmp = tmpl - coefl0(l)*pxl A 4 - coef 10 (2) *pxl A 3 
tmp = tmp - coef 11 (1) *pyl A 4 - coef 11 (2) *pyl A 3 - 
bl - [bl; tmp - coef 12 ( 1 ) *pzl A 4- coef 12 (2 ) *pzl A 



- coefl0(3) *pxl A 2 - coef 10 (4) *pxl A l; 

coefll (3) *pyl A 2 - coef 11 (4 ) *pyl A l; 
coefl2 (3) *pzl A 2 - coef 12 (4 ) *pzl A l; ] ; 



- coefl0(3) *pxl A 2 - coeflO (4) *pxl A l; 

coefll (3) *pyl A 2 - coef 11 ( 4 ) *pyl A l; 
3 - coef 12 (3) *pzl A 2 - coef 12 ( 4 ) *pzl A l; ] 



- coeflO (3) *pxl A 2 - 

- coefll (3) *pyl A 2 - 
3 - coefl2(3)*pzl A 2 



coeflO (4) *pxl A l; 
coefll (4) *pyl A l; 
- coefl2 (4) *pzl A l;] 



- coeflO (3) *pxl A 2 - coeflO (4) *pxl A l; 

coefll (3) *pyl A 2 - coef 11 (4 ) *pyl A l; 
3 - coefl2 (3) *pzl A 2 - coef 12 (4 ) *pzl A l; ] ; 



%; 
%; 
%; 



coefficients of the torus equation for (...) Tht = variable 



coef 13 = pinv(al)*bl; 

disp ( 1 coef 13 . . . 1 ) ; 

pxl = 2*rl; 

pyl = 2*rl; 

pzl = 2*rl; 

Tht = 0.0; 

Thtn = 3.14159/2; 

torusly = 0; 

for il = l:nr 

torusly = torusly + circlely; 
end 

al = [Thtn A 4 Thtn A 3 

tmp = tmpl - coeflO 
^:mp - tmp - coef 11 (1) *pyl / 
"mp = tmp - coef 12 (1) *pzl" 
"bl = [tmp - coef 13(1) *Tht' 

Thtn = 3.14159; 

torusly = 0; 

for il = l:nr 



Thtn A 2 Thtn A l;] ; 

l)*pxl A 4 - coeflO (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 ( 4 ) *pxl A l; 
4 - coefll (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4 ) *pyl A l; 
4- coefl2(2)*pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 (4) *pzl A l; 
4 - coefl3(2)*Tht A 3 - coef 13 (3) *Tht A 2 - coef 13 (4 ) *Tht A l; ] . 



torus 1^.- torusly + circlely; 
end 

al = [al; Thtn A 4 Thtn A 3 Thtn A 2 Thtn A l;]; 

tmp = tmpl - coeflO (1) *pxl A 4 : - - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 
"pp = tmp - coef 11 (1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4 ) *pyl A l; 
Imp = tmp - coefl2 (1) *pzl A 4- coef 12 (2) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 (4 ) *pzl A l; 
bl = [bl; tmp - coef 13 (1) *Tht A 4 - coef 13 (2) *Tht A 3 - coef 13 (3) *Tht A 2 - coef 13 (4 ) *Tht A l; ] ; 
Thtn = 3*3.14159/2; 
torusly = 0; 
for il = l:nr 

torusly - torusly + circlely; 
end 

al = [al; Thtn A 4 Thtn A 3 Thtn A 2 Thtn A l;]; 

tmp = tmpl - coeflO(l)*pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coeflO (4) *pxl A l; 
tmp = tmp - coefll (1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coefll (4) *pyl A l; 
tmp = tmp - coef 12 (1) *pzl A 4- coef 12 (2) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 ( 4 ) *pzl A l; 
bl = [bl; tmp - coef 13 (1) *Tht A 4 - coef 13 (2) *Tht A 3 - coef 13 (3) *Tht A 2 - coef 13 ( 4 ) *Tht A l; ] ; 
Thtn = 2*3.14159; 
torusly = 0; 
for il = l:nr 

torusly = torusly + circlely; 
end 

al = [al; Thtn A 4 Thtn A 3 Thtn A 2 Thtn A l;]; 

tmp = tmpl - coeflO (1) *pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 
tmp « tmp - coefll (1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 ( 3) *pyl A 2 - coef 11 ( 4 ) *pyl A l; 
tmp = tmp - coef 12 (1) *pzl A 4- coef 12 (2 ) *pzl A 3 - coef 12 ( 3) *pzl A 2 - coef 12 ( 4 ) *pzl A l ; 
bl = [bl; tmp - coef 13 (1) *Tht A 4 - coef 13 (2 ) *Tht A 3 - coef 13 (3) *Tht A 2 - coef 13 (4 ) *Tht A l; ] ; 

%; 

%; coefficients of the torus equation for (...) Thtn = variable 
%; 

coef 14 = pinv(al)*bl; "* vt 

iisp { 1 coef 14 . . . ? ) ; 

%; Now, finally the complete torus equation 
%; 

%;coef6 = eval(coef6); 
%;coef7 = eval(coef7); 
%;coef8 = eval(coef8); 
%;coef9 = eval (coef 9) ; 
%;coefl0 = eval (coeflO) 
%;coefll = eval (coefll) 
%;coefl2 = eval(coefl2) 
%;coefl3 = eval(coefl3) 
%;coef!4 = eval(coefl4) 

syms Rtt Otxt Otyt Otzt pxlt pylt pzlt thtt Thtnt; 

toruseqly = coef 6 ( 1) *Rtt A 4 + coef 6 (2 ) *Rtt A 3 + coef 6 (3) *Rtt A 2 + coef 6 (4 ) *Rtt A l; 

toruseqly = toruseqly + coef 7 ( 1) *Otxt A 4 + coef 7 (2) *Otxt A 3 + coef 7 (3) *Otxt A 2 + coef7(4) 
*Otxt A l; 

toruseqly = toruseqly + coef 8 (1) *Otyt A 4 + coef 8 (2) *Otyt A 3 + coef 8 (3) *Otyt A 2 + coef8(4) 
*Otyt A l; 

toruseqly = toruseqly + coef 9 ( 1 ) *Otzt A 4 + coef 9 (2 ) *Otzt A 3 + coef 9 (3) *Otzt A 2 + coef 9 (4) 
*0tzt A l; 

toruseqly = toruseqly + coef 10 (1) *pxlt A 4 + coef 10 (2) *pxlt A 3 + coef 10 (3) *pxlt A 2 + coeflO (4) 
*pxlt A l; 

toruseqly = toruseqly + coef 11 ( 1 ) *pylt A 4 + coef 11 (2) *pylt A 3 + coef 11 (3) *pylt A 2 + coefll (4) 
!?pylt A l; 

Joruseqly = toruseqly + coef 12 ( 1 ) *pzlt A 4 + coef 12 (2) *pzlt A 3 + coef 12 (3) *pzlt A 2 + coef 12 (4) 
p *pzlt A l; 

toruseqly = toruseqly + coef 13 (1) *Thtt A 4 + coef 13 (2) *Thtt A 3 + coef 13 (3) *Thtt A 2 + coef 13 (4) 
*Thtt A l; 

toruseqly - toruseqly.- -Kcoef 14 (1) *Thtnt A 4 + coef 14 (2 ) *ThtnftS3 + coef 14 (3) *Thtnt A 2 + coef 14 



(4)*Thtnt A l; 

disp ( 1 finished computation of toruseqly 
fprintf (f id, 1 \n\ntoruseqly = Sum of: \n 
for il - 1:4 

fprintf (fid, ' %20 . 12f ' , coef 6 (il) ) ; 

fprintf (fid, 1 Rtt A %li ' , il) ; 

f print f (fid, '\n' ) ; 
end 

for il = 1:4 

fprintf (f id, • %20 . 12f 1 , coef 7 (il) ) ; 

fprintf (fid, ' Otxt A %li 1 , il) ; 

fprintf (fid, 1 \n' ) ; 
end 

for il = 1:4 

fprintf (fid, 1 %20 . 12f ' , coef 8 ( il) ) ; 

fprintf (fid, 1 Otyt A %li 1 , il) ; 

fprintf (fid, '\n' ) ; 
end 

for il = 1:4 

fprintf (fid, ' %20 . 12f 1 , coef 9 (il ) ) ; 

fprintf (fid, 1 Otzt A %li ' , il) ; 

fprintf (fid, f \n' ) ; 
end 

for il = 1:4 

fprintf (fid, ' %20 . 12f 1 , coef 10 (il) ) ; 

fprintf (fid, 1 pxlt A %li f , il) ; 

fprintf (fid, f \n' ) ; 
end 

for il = 1:4 

fprintf (fid, 1 %20 . 12f ? , coef 11 (il) ) ; 
fprintf (fid, ' pylt A %li ■ , iT) ; 
fprintf (fid, '\n f ) ; 

Pnd 
or il = 1:4 
fprintf (fid, 1 %20 . 12f 1 , coef 12 (il) ) ; 
fprintf (fid, ' pzlt A %li 1 , il ) ; 
fprintf (fid, '\n f ) ; " 
end 

for il - 1:4 

fprintrfTKd-, r: %20 . 12f 1 , coefl3 (il) ) ; 

fprintf (fid, ' Thtt A %li 1 , il) ; 

fprintf (fid, '\n ? ) ; 
end 

for il = 1:4 

fprintf (fid, » %20 . 12f 1 , coef 14 (il) ) ; 

fprintf (fid, 1 Thtnt A %li 1 , il ) ; 

fprintf (fid, 1 \n' ) ; 
end 



fprintf (fid, 1 \n\nTest of torusly 1 ); 

Rtt = 5.0; 

Otxt = 0.0; 

Otyt = 0.0; 

Otzt = 0.0; 

pxlt = 0.0; 

pylt = 0.0; 

pzlt = 0.0; 

Thtt = 0.0; 

Thtnt = 0.0; 

fprintf (fid, 1 \nRtt = %20 • 12f ' , Rtt ) ; 

I fprintf (fid, ' \nOtx = %20 . 12f 1 , Otxt ) ; 
fprintf (fid, 1 \nOty = %20 . 12f ' , Otyt ) ; 
fprintf (fid, f \nOtz « %20. 12f 1 , Otzt) ; 
fprintf (fid, 1 \npxlt = %20. 12f 1 ,pxlt) ; 
fprintf (fid, f \npylt = %2£lU2f ' , pylt ) ; 
fprintf (fid, 1 \npzlt = %20~ T2f 1 , pzlt ) ; 



f print f (fid, 
f print f (fid, 
f print f (fid, 
fprintf (f rdr, 
'fprintf (fid, 
\tt = 5.0; 



'\nThtt = %20.12f ',Thtt) ; 
•\nThtnt = %20 . 12 f 1 , Thtnt ) ; 
'\ntorusly = %20 . 12f 1 , eval (torusly) ) 

1 \n\nTest of torusly 1 ) ; 



Otxt = 
Otyt » 
Otzt = 
pxlt = 
pylt - 
pzlt = 
Thtt = 
Thtnt = 



0.0; 
0.6; 
0.0; 
10.0; 
0.0 
0.0 
0.0 
0.0; 



fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
Rtt = 5.0; 
Otxt =0.0 
Otyt =0.0 
Otzt = 0.0; 
pxlt = 0.&; 
pylt =0.0; 
)zlt = 10.0; 
f-httr = 
"Thtnt = 0.0; 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 



' \nRtt = %20.12f ' ,Rtt) ; 

r \nOtx = %20.12f *,Otxt) ; 

f \nOty = %20.12f 1 ,Otyt) ; 

'\nOtz = %20.12f 1 ,Otzt) ; 

'\npxlt = %20.12f ',pxlt) ; 

'\npylt = %20.12f 1 ,pylt) ; 

'\npzlt = %20.12f 1 ,pzlt) ; 

'\nThtt = %20.12f 1 f Thtt) ; 

'\nThtnt = %20. 12f ! , Thtnt) ; 

'\ntorusly = %20 . 12f 1 , eval (torusly) ) ; 

r \n*); 

r \n\nTest of torusly 1 ); 



\nRtt = %20, 
\nOtx = %20. 
\nOty = %20. 
\nOtz = %20. 



12f f ,Rtt) ; 
12f 1 ,Otxt) , 
12f',0tyt) 
12f ? ,Otzt) 



Rtt 
Otxt 
Otyt 
Otzt 
pxlt 
pylt 
pzlt 
Thtt 



= 5.0; 
0.0 
0.0 
0.0 
-10 
0.0 
0.0 
0.0 



\npxlt = %20.12f 1 ,pxlt) ; 
\npylt = %20.12f 1 ,pylt) ; 
\npzlt = %20.12f 1 ,pzlt) ; 
\nThtt = %20.12f 1 ,Thtt) ; 
\nThtnt = %20 . 12f 1 , Thtnt ) ; 
\ntorusly = %20 . 12f 1 , eval (torusly) ) 
An' ) ; 

\n\nTest of torusly 1 ); 



Thtnt = 0.0; 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 

^fprintf (fid, 
fprintf (fid, 

^fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 



\nRtt 
\nOtx 
\nOty 
\nOtz 
\npxlt 
\npylt 
\npzlt 
\nThtt 
\nThtnt 



= %20, 
= %20. 
= %20. 
= %20. 



12f 1 ,Rtt) ; 
12f 1 ,Otxt) ; 
12f',Otyt); 
12f 1 ,Otzt) ; 



= %20.12f f ,pxlt) ; 
= %20.12f 1 ,pylt) ; 
= %20.12f f ,pzlt) ; 
= %20.12f f ,Thtt); 
= %2JD. 12 f», Thtnt) 



\ntorusly = %20 . 12f 1 , eval (torusly) ) 



• 



fprintf (fid, '\n') ; 

fprintf (fid, f \n\nTest of torusly'); 
Rtt = 5.0; 
Otxt = Org?" 
fcyt = OTO; 
tzt =0.0; ~ 
pxlt = 0.0; 
pylt = 0.0; 
pzlt = -10.0; 
Thtt = 0.0"; 
Thtnt = 0.0; 

fprintf (fid, 1 \nRtt « %20. 12f 1 , Rtt) ; 

fprintf (fid, 1 \nOtx = %20 . 12f f , Otxt ) ; 

fprintf (fid, f \nOty = %20 . 12f 1 , Otyt ) ; 

fprintf (fid, f \nOtz = %20. 12f 1 ,Otzt) ; 

fprintf (fid, f \npxlt = %20 . 12f \ pxlt ) ; 

fprintf (fid, 1 \npylt = %20 . 12f ' , pylt ) ; 

fprintf (fid, '\npzlt = %20 . 12f • , pzlt ) ; 

fprintf (fid, f \nThtt = %20.12f 1 ,Thtt) ; 

fprintf (fid, '\nThtnt = %20 . 12f f , Thtnt ) ; 

fprintf (fid, f \ntorusly = %20 . 12f 1 , eval (torusly) ) ; 

fprintf (fid, f \n f ); 

syms coef6 coef7 coef8 coef9 coeflO coefll coefl2 coefl3 coefl4; 

rl = 1.0; 
Otx = 0; 
Oty = 0; 
Otz = 0; 

xpt = Otx + Rt*cos(Tht); 
ypt = Oty + Rt*sin(Tht) ; 

rpt = sqrt((xpt - Otx) A 2 + (ypt - Oty) A 2); 
t = xpt; 
ypt; 

%; thpt = itanl (xt , yt ) ; 
syms thpt; 

thcum = thpt + Thtn; 
xO = rpt^cos (thcum) ; 
yO = rpt*sin (thcum) ; 
thl = Tht; 
th2 = 0; 
nr = 8; 

thp = 2*3.14159/nr; 
Tht = thp*JLl; 
Rt = 20; 
toruslz = 0; 
for il » l:nr 

thpt = itanl (xt, yt) ; 

toruslz = toruslz + circlelz; 
end 

al « [ Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 
bl = [eval (toruslz) ;] ; 
Rt = 60; 
toruslz = 0; 
for il = l:nr 

thpt = itanl (xt, yt) ; 

toruslz = toruslz + circlelz; 
end 

al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 
bl = [bl; eval (toruslz) ;] ; 
t - 240; 
oruslz = 0; 
'for il = l:nr 

thpt = itanl (xt, yt) ; 
toruslz « toruslz + circlelz; 
end 



al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 2;]; 
bl = [bl; eval (toruslz) ; ] ; 
Rt = 2400; 
toruslz = 0; 
or il = l:nr 
thpt = itanl (xt , yt ) ; 
toruslz;= trortrslz + circlelz; 
end 

al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 2 ; ] ; 
bl = [bl; eval (toruslz) ;] ; 
% ; 

%; The following is the first set of coefficients of the torus equation for Rt = variable 
ft; 

coef6 = pinv(al)*bl; 

disp( 1 coef 6. . . ? ) ; 

xO = Otx + Rt*cos(Tht); 

yO = Oty + Rt*sin(Tht); 

thl = Tht; 

nr = 8; 

thp = 2*3.14159/nr; 
Tht = thp*il; 
Rt = 20; 
Otx = 8*rl; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [Otx A 4 Otx A 3 Otx A 2 Otx A l;]; 

bl = [ (eval (toruslz) - coef 6 ( 1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l) ; ] ; 
Otx = 64*rl; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
nd 

al = [al; Otx A 4 Otx A 3 Otx A 2 Otx A l;]; 

bl = [bl; (eval (toruslz) - coef 6 (1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef6(4)*Rt A 
1);]; 

Otx = -8*rl; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [al; Otx A 4 Otx A 3 Otx A 2 Otx A l;]; 

bl = [bl; (eval (toruslz) - coef 6 (1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef6(4)*Rt A 
1);]; 

Otx = -64*rl; 
toruslz = 0; 
for il = l:nr 

toruslz - toruslz + circlelz; 
end 

al - [al; Otx A 4 Otx A 3 Otx A 2 Otx A l;]; 

bl = [bl; (eval (toruslz) - coef 6 (1) *Rt A 4 - coef 6 (2 ) *Rt A 3 - coef 6 (3) *Rt A 2 - coef6(4)*Rt A 
l);]; 

ft; 

%; The following is the first set of coefficients of the torus equation for Rt = 20, Otx 

variable 

ft; 

coef 7 = pinv(al)*bl; 
disp( 1 coef 7 ...'); 
0 = Otx + Rt*cos (Tht) ; 
0 = Oty + Rt*sin(Tht) ; 
hi = Tht; 
nr = 8; 

thp = 2*3.14159/nr; 

Tht = thp*il; ..:::...;•» 



Rt = 20^: ; 
Otx = 0; 
Oty = 8*rl; 
toruslz =^0; 

•or il = l:nr 
toruslz = toruslz + circlelz; 
end 

al = [Oty A 4 Oty A 3 Oty A 2 Oty A l; ] ; 

tmp = eval (toruslz) - coef 6 (1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l; 
bl = [tmp - coef7 (l)*Otx A 4 - coef 7 (2 ) *Otx A 3 - coef7 (3) *Otx A 2 - coef 7 (4 ) *0tx A l; ] ; 
Oty = 64*rl; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [airOty A 4 Oty A 3 Oty A 2 Oty A l;]; 

tmp = eval (toruslz) - coef 6 (1) *Rt A 4 - coef 6 (2 ) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4) *Rt A l; 
bl = [bl; tmp - coef 7 (1) *Otx A 4 - coef 7 (2) *Otx A 3 - coef 7 (3) *0tx A 2 - coef7 (4) *Otx A l; ] ; 
Oty = -8*rl; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [al; Oty A 4 Oty A 3 Oty A 2 Oty A l;]; 

tmp = eval (toruslz) - coef 6 (1) *Rt A 4 - coef 6 (2 ) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l; 
bl = [bl; tmp - coef 7 (1) *Otx A 4 - coef 7 (2) *Otx A 3 - coef 7 (3) *0tx A 2 - coef7 (4) *0tx A l; ] ; 
Oty = -64*rl; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al - [al; 0ty A 4 Oty A 3 Oty A 2 Oty A l;]; 

•mp = eval (toruslz) - coef 6 ( 1) *Rt A 4 - coef 6 (2 ) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l; 
1 - [bl; tmp - coef7 (1) *Otx A 4 - coef 7 (2 ) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 (4 ) *Otx A l; ] ; 

%; 

%; The following is the first set of coefficients of the torus equation for Rt = 20, Otx 

0, Oty = variable 

%; 

coef8 = pinv(al)*bl; 

disp( 'coef 8. . • 1 ) ; 

xO = Otx + Rt*cos(Tht); 

yO = Oty + Rt*sin(Tht); 

thl = Tht; 

nr = 8; 

thp = 2*3.14159/nr; 
Tht = thp*il; 
Rt = 20; 
Otx = 0; 
Oty = 0; 
Otz = 8*rl; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [Otz A 4 Otz A 3 Otz A 2 Otz A l;]; 

tmp = eval (toruslz) - coef 6 (1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l; 
tmp = tmp - coef7(l)*Otx A 4 - coef 7 (2) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 (4) *Otx A l; 
bl = [tmp - coef8(l)*Oty A 4 - coef 8 (2 ) *Oty A 3 - coef 8 (3) *Oty A 2 - coef 8 (4 ) *Oty A l; ] ; 
Otz = 64*rl; 

•toruslz = 0; 
for il = l:nr 
toruslz = toruslz + circlelz; 
end 

al = [al; Otz A 4 Otz A 3 Otz A 2 Otz A l;]; 

tmp = eval (toruslz) - coef 6 { 1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 ( 4 ) *Rt A l; 



tmp = tmp - coef7 (1) *0tx A 4 - coef 7 (2) *0tx A 3 - coef 7 (3) *0tx A 2 - coef7 (4) *Otx A l; 

bl = [bl; tmp - coef 8 (1) *0ty A 4 - coef 8 (2) *0ty A 3 - coef 8 (3) *0ty A 2 - coef 8 (4 ) *Oty A l; ] ; 

Otz = -8*rl; 

toruslz = 0; 

•or il = l:nr 
toruslz = toruslz + circlelz; 
end 

al = [al; 0tz A 4 Otz A 3 Otz A 2 Otz A l;]; 

tmp = eval (toruslz) - coef 6 (1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l; 

tmp = tmp - coef7 (1) *Otx A 4 - coef 7 (2) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 (4) *Otx A l; 

bl = [bl; tmp - coef 8 ( 1 ) *Oty A 4 - coef 8 (2) *Oty A 3 - coef 8 (3) *Oty A 2 - coef 8 ( 4 ) *Oty A l; ] ; 

Otz = -64*rl; 

toruslz = 0; 

for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [al; Otz A 4 Otz A 3 Otz A 2 Otz A l;]; 

tmp = eval (toruslz) - coef 6 (1) *Rt A 4 - coef 6 (2 ) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l; 

tmp = tmp - coef7 (1) *Otx A 4 - coef 7 (2) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 (4 ) *Otx A l; 

bl = [bl; tmp - coef 8 (1) *Oty A 4 - coef 8 (2) *Oty A 3 - coef 8 (3) *Oty A 2 - coef 8 (4 ) *0ty A l; ] ; 

%; 

%; The following is the first set of coefficients of the torus equation for Rt = 20, Otx = 

0, Oty = 0, Otz = var. 

%; 

coef 9 = pinv(al)*bl; 
disp{ 1 coef 9. . . 1 ) ; 
%; 

xO = Otx + Rt*cos(Tht); 
yO = Oty + Rt*sin(Tht); 
thl = Tht; 
nr = 8; 

•hp = 2*3.14159/nr; 
ht - thp*il; 
Rt = kl*rl; 
Otx = k2*rl; 
Oty = k3*rl; 
Otz = k4*rl; 
pxl = 2*rl; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 

tmpl = eval (toruslz) - coef 6 (1) *Rt A 4 - coef 6 (2) *Rt A 3 - coef 6 (3) *Rt A 2 - coef 6 (4 ) *Rt A l; 
tmpl = tmpl - coef7 (l)*Otx A 4 - coef 7 (2) *Otx A 3 - coef 7 (3) *Otx A 2 - coef 7 ( 4 ) *Otx A l; 
tmpl = tmpl - coef8 (1) *Oty A 4 - coef 8 (2 ) *Oty A 3 - coef 8 (3) *Oty A 2 - coef 8 ( 4 ) *Oty A l; 
tmpl = tmpl ~ coef9(l)*Otz A 4 - coef 9 (2) *Otz A 3 - coef 9 (3) *Otz A 2 - coef 9 (4 ) *Otz A l; 



bl = [tmpl; ] ; 
pxl = 8*rl; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [al; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 
bl = [bl; tmpl; ] ; 
pxl = 32*rl; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
^nd 

"al = [al; pxl A 4 pxl A 3 pxl A 2 pxl A l;]; 
bl = [bl; tmpl;] ; 
pxl = 200*rl; 
toruslz = 0; 



t 



for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [ales**** pxl A 3 pxl A 2 r; pxl A l; ] ; 
"|1 = [bl; tmpl;]; 

%; coefficients of the torus equation for (...) pxl = variable 
%; 

coeflO = pinv(al)*bl; 
disp( 'coeflO. . . 1 ) ; 

pxl = 2*rl; 
pyl = 2*rl; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

bl = [tmpl - coeflO (1) *pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4) *pxl A l; ] / 
pyl = 8*rl; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [al; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

bl = [bl; tmpl - coef 10 (1) *pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 ( 3) *pxl A 2 - coef 10 (4) *pxl A l; ] ; 
pyl = 32*rl; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [al; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 
1 = [bl; tmpl - coeflO (1) *pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; ] ; 
yl = 200*rl; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [al; pyl A 4 pyl A 3 pyl A 2 pyl A l;]; 

bl = [bl; tmpl - coef 10 (1) *pxl A 4 - coef 10 (2 ) *pxl A 3 - coef 10 (3) *pxl A 2 - coeflO (4) *pxl A l; ] ; 
%; 

coefficients of the torus equation for (...) pyl = variable 



coefll = pinv(al)*bl; 
disp( 'coefll. . . 1 ) ; 
pxl = 2*rl; 
pyl = 2*rl; 
pzl = 2*rl; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

tmp = tmpl - coeflO (l)*pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 
bl = [tmp - coefll (l)*pyl A 4 - coef 11 (2 ) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4 ) *pyl A l; ] ; 
pzl = 8*rl; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

k al = [al; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

:mp = tmpl - coef 10 (1) *pxl A 4 - coef 10 (2 ) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 
'bl - [bl; tmp - coef 11 (1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4) *pyl A l; ] ; 

pzl = 32*rl; 

toruslz = 0; 

for il = l:nr 



toruslz = toruslz + circlelz; 
end 

al = [al; pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 



tmp = tmpl - coeflO(l)*pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 

«1 = [bl; tmp - coefll (l)*pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coefll (4) *pyl A l; ] ; 
zl = 200*rl; 
toruslz = -0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [al; "pzl A 4 pzl A 3 pzl A 2 pzl A l;]; 

tmp = tmpl - coefl0(l)*pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 

bl = [bl; tmp - coef 11 (1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 ( 4 ) *pyl A l; ] ; 

%; 

%; coefficients of the torus equation for (...) pzl = variable 
%; 

coef 12 = pinv(al)*bl; 

disp( 1 coef 12. . . 1 ) ; 

pxl = 2*rl; 

pyl = 2*rl; 

pzl = 2*rl; 

Tht = 3.14159/2; 

toruslz = 0; 

for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [Tht A 4 Tht A 3 Tht A 2 Tht A l;]; 

tmp = tmpl - coefl0(l)*pxl A 4 - coef 10 (2 ) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 
tmp = tmp - coefll (1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4 ) *pyl A l; 
bl = [tmp - coe*12(l)*pzl A f- coef 12 (2) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 (4 ) *pzl A l; ] ; 
Tht = 3.14159; 
^toruslz == 0; 
tor il = l-:nr 

toruslz = toruslz + circlelz; 
end 

al = [al; Tht A 4 Tht A 3 Tht A r Tht A l; ] ; 

tmp = tmpl - coeflO(l)*pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 
tmp = tmp - coefll (l)*pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3 ) *pyl A 2 - coef 11 ( 4 ) *pyl A l ; 
bl = [bl; tmp - coef 12 (1) *pzl A 4- coef 12 (2) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 (4 ) *pzl A l; ] ; 
Tht = 3*3.14159/2; 
toruslz == 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [al; Tht A 4 Tht A 3 Tht A 2 Tht A l;]; 

tmp = tmpl - coefl0(l)*pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 
tmp = tmp - coefll (1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4 ) *pyl A l; 
bl = [bl; tmp - coef 12 (1) *pzl A 4- coef 12 (2 ) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 (4 ) *pzl A l; ] ; 
Tht = 2*3.14159; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [al; Tht A 4 Tht A 3 Tht A 2 Tht A l;]; 

tmp - tmpl - coef!0(l)*pxl A 4 - coef 10 (2 ) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 
tmp = tmp - coefll (1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4) *pyl A l; 
bl « [bl; tmp - coef 12 ( 1) *pzl A 4- coef 12 (2 ) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 (4 ) *pzl A l; ] ; 



m 



coefficients of the torus equation for (...) Tht = variable 



coefl3 = pinv(al)*bl; 
disp( 'coef 13. . . ? ) ; 
pxl = 2*rl; 
pyl = 2*rl; 



pzl = 2*rl; *' " s "~- 

Tht = 0.6;' 
Thtn = 3.14159/2; 
torus 1 z =-■ 07 — 
Aor il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [Thtn A 4 Thtn A 3 Thtn A 2 Thtn A l;]; 

tmp = tmpl - coefl0(l)*pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4) *pxl A l; 
tmp - tmp - coefll(l)*pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4 ) *pyl A l; 
tmp = tmp - coefl2 (1) *pzl A 4- coef 12 (2) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 (4) *pzl A l; 
bl = [tmp - coefl3(l)*Tht A 4 - coef 13 (2) *Tht A 3 - coef 13 (3) *Tht A 2 - coef 13 (4 ) *Tht A l; ] ; 
Thtn = 3.14159; 
toruslz = 0; 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [al; Thtn A 4 Thtn A 3 Thtn A 2 Thtn A l;]; 

tmp = tmpl - coeflO (1) *pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 

tmp = tmp - coefll (1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4 ) *pyl A l; 

tmp = tmp - coef 12 (1) *pzl A 4- coef 12 (2) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 (4) *pzl A l; 

bl = [bl; tmp - coef 13 (1 ) *Tht A 4 - coef 13 (2) *Tht A 3 - coef 13 (3) *Tht A 2 - coef 13 (4 ) *Tht A l; ] ; 

Thtn = 3*3.14159/2; 

toruslz = 0; 

for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [air^l^rtTT^ Thtn A 3 Thtir^ Thtn A l;]; 

tmp = tmpl - coef 10(1) *pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 ( 3) *pxl A 2 - coef 10 (4) *pxl A l; 
tmp = tmp - coefll (1) *pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4 ) *pyl A l; 
tmp = tmp - coefl2 (1) *pzl A 4- coef 12 (2 ) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 (4 ) *pzl A l; 
bl = [bl; tmp - coef 13 (1) *Tht A 4 - coef 13 (2 ) *Tht A 3 - coef 13 (3) *Tht A 2 - coef 13 (4) *Tht A l; ] ; 

•htn = 2*3.14159; 
oruslz =-0r 
for il = l:nr 

toruslz = toruslz + circlelz; 
end 

al = [al; Thtn A 4 Thtn A 3 Thtn A 2 Thtn A l;]; 

tmp - tmpl - coef 10 (1) *pxl A 4 - coef 10 (2) *pxl A 3 - coef 10 (3) *pxl A 2 - coef 10 (4 ) *pxl A l; 

tmp = tmp - coefll (l)*pyl A 4 - coef 11 (2) *pyl A 3 - coef 11 (3) *pyl A 2 - coef 11 (4 ) *pyl A l; 

tmp = tmp - coefl2(l)*pzl A 4- coef 12 (2 ) *pzl A 3 - coef 12 (3) *pzl A 2 - coef 12 (4) *pzl A l; 

bl = [bl; tmp - coef 13 (1) *Tht A 4 - coef 13 (2 ) *Tht A 3 - coef 13 (3) *Tht A 2 - coef 13 (4 ) *Tht A l; ] ; 



%; 

%; coefficients of the torus equation for (...) Thtn = variable 
%; 

coefl4 = pinv(al)*bl; 
disp ( 1 coef 14 . . . 1 ) ; 

%; 

%; Now, finally the complete torus equation 
%; 

%;coef6 = eval(coef6); 
%;coef7 = eval(coef7); 
%;coef8 = eval(coef8); 
%;coef9 = eval(coef9); 
%; coef 10 = eval(coeflO) 
%;coefll = eval (coefll) 
%;coefl2 = eval (coef 12) 

I%;coefl3 = eval (coef 13) 
|;coefl4 = eval(coefl4) 

syms Rtt Otxt Otyt Otzt pxlt pylt pzlt Thtt Thtnt; 



toruseqlz = coef 6 ( 1 ) *Rtt A 4 + coef 6 (2 ) *Rtt A 3 + coef 6 (3) *Rtt A 2 + coef 6 ( 4 ) *Rtt A l ; 



toruseqlz = toruseqlz + coef7 (1) *0txt A 4 + coef 7 (2) *0txt A 3 + coef 7 (3) *0txt A 2 + coef7(4) 
*Otxt A l; 

toruseqlz = toruseqlz + coef8 (1) *0tyt A 4 + coef 8 (2) *Otyt A 3 + coef 8 (3) *0tyt A 2 + coef8(4) 
*Otyt A l; 

oruseqlz = toruseqlz + coef 9 ( 1 ) *0tzt A 4 + coef 9 (2) *Otzt A 3 + coef 9 ( 3) *Otzt A 2 + coef9(4) 
Otzt A l; 

toruseqlz = toruseqlz + coef 10 ( 1 ) *pxlt A 4 + coef 10 (2) *pxlt A 3 + coef 10 (3) *pxlt A 2 + coef 10 (4) 
*pxlt A l; 

toruseqlz = toruseqlz + coef 11 ( 1) *pylt A 4 + coef 11 (2 ) *pylt A 3 + coef 11 (3) *pylt A 2 + coefll{4) 
*pylt A l; . 

toruseqlz = toruseqlz + coef 12 ( 1) *pzlt A 4 + coef 12 (2 ) *pzlt A 3 + coef 12 ( 3) *pzlt A 2 + coef 12 (4) 
*pzlt A l; 

toruseqlz = toruseqlz + coef 13 ( 1) *Thtt A 4 + coef 13 (2 ) *Thtt A 3 + coef 13 (3) *Thtt A 2 + coef 13 (4) 
*Thtt A l; 

toruseqlz = toruseqlz + coef 14 ( 1) *Thtnt A 4 + coef 14 (2) *Thtnt A 3 + coef 14 (3) *Thtnt A 2 + coef!4 
(4)*Thtnt A l; 

disp ( 1 finished computation of toruseqlz . • . 1 ) ; 
fprintf (fid, 1 \n\ntoruseqlz = Sum of: \n f ); 
for il = 1:4 

fprintf (fid, ' %20 . 12f ' , coef 6 (il) ) ; 

fprintf (fid, 1 Rtt A %li f , il) ; 

fprintf (fid, 'Xn 1 ) ; 
end 

for il = 1:4 

fprintf (fid, 1 %20. 12f f , coef7 (il) ) ; 

fprintf (fid, 1 Otxt A %li f , il) ; 

fprintf (fid, f \n f ) ; 
end 

for il = 1:4 

fprintf ( fid, 1 %20 . 12f f , coef 8 ( il ) ) ; 
fprintf (fid, 1 Otyt A %li 1 , il) ; 
fprintf (fid, f \n' ) ; 
^end 
for il = 1:4 

fprintf (fid, 1 %20 . 12f 1 , coef 9 ( il ) ) ; 
fprintf ( fid, Otzt A %li 1 , il) ; 
fprintf (fid, f \n' ) ; 
end 

for il « 1:4 

fprintf (fid, 1 %20 . 12f f , coef 10 (il ) ) ; 

fprintf (fid, ' pxlt A %li 1 , il) ; 

fprintf (fid, , \ri r ) ; 
end ^ 
for il = 1:4 

fprintf (fid, 1 %20.12f 1 , coef 11 (il) ) ; 

fprintf (fid, r pylt A %li 1 , il ) ; 

fprintf (fid, f \n ? ) ; 
end 

for il = 1:4 

fprintf (fid, 1 %20.12f 1 , coef 12 (il) ) ; 

fprintf (fid, f pzlt A %li f , il) ; 

fprintf (fid, 1 \n» ) ; 
end 

for il = 1:4 

fprintf (fid, 1 %20 . 12f f , coef 13 (il) ) ; 

fprintf (fid, 1 Thtt A %li 1 , il) ; 

fprintf (fid, 'Xn') ; 
end 

for il = 1:4 

fprintf (fid, ? %20. 12f 1 , coefl4 (il) ) ; 
fprintf (fid, 1 Thtnt A %li 1 , il) ; 
fprintf (fid, 1 \n ! ) ; 
r end 

fprintf (fid, ' \n\nTest of toruslz 1 ); 
Rtt =5.0; 



Otxt = 0.0; 

Otyt = 0.0; 

Otzt = 0.0; 

pxlt = O/gf 
^fcpylt = 0.6; 
^m>zlt = 0.0; 

Thtt = 0.0; 

Thtnt = 0.0; 

fprintf (fid, 1 \nRtt = %20 . 12f f , Rtt ) ; 

fprintf (fid, 1 \nOtx = %20 . 12f 1 , Otxt ) ; 

fprintf (fid, 1 \nOty = %20 . 12f ' , Otyt ) ; 

fprintf (fid, '\nOtz = %20 . 12f 1 , Otzt ) ; 

fprintf (fid, 1 \npxlt = %20. 12f 1 , pxlt) ; 

fprintf (fid, 1 \npylt = %20 . 12f f , pylt ) ; 

fprintf (fid, '\npzlt - %20 . 12f ' , pzlt ) ; 

fprintf (fid, 1 \nThtt = %20 . 12f 1 , Thtt ) ; 

fprintf (fid, 1 \nThtnt = %20 . 12f , Thtnt ) ; 

fprintf (fid, f \ntoruslz = %20 . 12f , eval (toruslz) ) ; 

fprintf (fid, »\n f ) ; 

fprintf (fid, 1 \n\nTest of toruslz 1 ); 



Rtt = 


5.0; 


Otxt = 


0.0; 


Otyt = 


0.0; 


Otzt = 


0.0; 


pxlt = 


10.0; 


pylt = 


0.0; 


pzlt = 


0.0; 


Thtt = 


0.0; 


Thtnt 


= 0.0; 



fprintf (fid, 1 \nRtt = %20.12f 1 ,Rtt) ; 
fprintf (fxd, f \nOtx = %20.12T f ,Otxt) ; 
fprintf (fid, 1 \nOty = %20 . 12f f , Otyt ) ; 

I fprintf (fid, 1 \nOtz = %20.12f 1 ,Otzt) ; 
fcprintf {f*a, 1 Wxlt =^20."S2f 1 ,pxlt) ; 
fprintf (fid, ' \npylt = %20 . 12f 1 , pylt ) ; 
fprintf (fid, '\npzlt = %20. 12f 1 ,pzlt) ; 
fprintf (fid, 1 \nThtt = %20 . 12f 1 , Thtt ) ; 
fprintf (fid, 1 \nThtnt = %20 . 12f 1 , Thtnt ) ; 
fprintf (fid, 'Xntoruslz = %20 . 12f , eval (toruslz) ) ; 
fprintf (fid, ? \n f ) ; 

fprintf (fid, 1 \n\nTest of toruslz 1 ); 

Rtt = 5.0; 

Otxt = 0.0; 

Otyt = 0.0; 

Otzt = 0.0; 

pxlt = 0.0; 

pylt = 0.0; 

pzlt = 10.0; 

Thtt = 0.0; 

Thtnt = 0.0; 

fprintf (fid, 1 \nRtt = %20.12f 1 ,Rtt) ; 

fprintf (fid, f \nOtx = %20.12f 1 ,Otxt) ; 

fprintf (fid, 1 \nOty = %20 . 12f 1 , Otyt ) ; 

fprintf (fid, f \nOtz = %20.12f ? ,Otzt) ; 

fprintf (fid, f \npxlt = %20 . 12f ? , pxlt ) ; 

fprintf (fid, '\npylt = %20. 12f 1 , pylt) ; 

fprintf (fid, '\npzlt = %20.12f 1 ,pzlt) ; 

fprintf (fid, ! \nThtt = %20 . 12f 1 , Thtt ) ; 

fprintf (fid, f \nThtnt = %20.12f 1 , Thtnt) ; 

fprintf (fid, '\ntoruslz = %20. 12f f , eval (toruslz) ) ; 

fprintf (fid, ? \n' ) ; 

^fprintf (fid, 1 \n\nTest of toruslz 1 ); 
Rtt = 5.0; 
Otxt = 0.0; 
Otyt = 0.0; 
Otzt = 0.0; 



in 

m 



pxlt = -10 
pylt =0.0 
pzlt = 0.0 
Thtt =0.0 

tnt = 0.0; 
'print f (fid, 
fprintf (f id, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
Rtt = 5.0; 
Otxt =0.0 
Otyt =0.0 
Otzt = 
pxlt = 
pylt = 
pzlt = 
Thtt = 
Thtnt = 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 
print f (fid, 
^ print f (fid, 
fprintf (fid, 
fprintf (fid, 
fprintf (fid, 



0; 



f \nRtt = %20.12f \Rtt) ; 

1 \nOtx = %20.12f ' ,Otxt) ; 

'\nOty = %20.12f f ,Otyt) ; 

'\nOtz = %20.12f \Otzt) ; 

'\npxlt = %20.12f ' ,pxlt) ; 

'\npylt = %20.12f \pylt) ; 

'\npzlt = %20.12f f ,pzlt) ; 

'\nThtt = %20.12f \Thtt); 

' \nThtnt = %20.12f Thtnt ) ; 

'\ntoruslz = %20 . 12f ' , eval (toruslz) ) ; 

'\n»); 

'\n\nTest of toruslz 1 ); 



0.0 
0.0 
0.0 
-10 
0.0 
0 



0; 



0; 



'\nRtt * %20.12f ? ,Rtt) ; 

'\nOtx = %20.12f 1 ,Otxt) ; 

'\nOty = %20.12f »,Otyt) ; 

'\nOtz = %20.12f 1 ,Otzt) ; 

'Vhpxlt = %20.1"2f 1 ,pxlt) ; 

( \npylt = %20.12f 1 ,pylt) ; 

'\npzlt = %20.12f f ,pzlt) ; 

1 \nThtt = %20.12f f ,Thtt) ; 

'\nThtnt = %20 . 12 f 1 , Thtnt ) ; 

1 \ntoruslz = %20 . 12f f , eval (toruslz) ) ; 

'\n'); 



compute the multi-torus equation 

- . . . f (x,y, z,rl,r2, thl, th2,r3,thml) 

- cumulative equation generated with a for loop 



pxlt 
pylt 
pzlt 
thlt 
th2t 



pxl; 

pyi; 

pzl; 
thl; 
th2; 



syms Rmt Thmt nT tmpmTx tmpmTy tmpmTz coeflS coefl6; 

Otxt = Rmt *cos (Thmt) ; 

Otyt = Rmt*sin (Thmt) ; 

Otzt = 0.0; 

Thtt = 0.0; 

Thtnt = 0.0; 

nT = 16; 

Rtt = 20*rl; 

•tmpmTx = 0; 
■or il = l:nT 
Thmt = il*2*3.14159/nT; 
tmpmTx = tmpmTx + toruseqlx; 
end 



Rmt = 


100Q*rl; 




al = 


[Rmt A 4 Rmt A 3 Rmt A 2 Rmt' 


s l; ] ; 


bl = 


[eval (tmpmTx) ; ] ; 
200.Q*rlr 




Rmt = 




1? = 


[al; Rmt A 4 Rmt A 3 Rmt A 2 


Rmt A l; ] ; 


ml = 


[bl; eval (tmpmTx) ; ] ; 




Rmt = 


10000*rl; 




al = 


[al; Rmt A 4 Rmt A 3 Rmt A 2 


Rmt A l; ] ; 


bl = 


[bl; eval (tmpmTx) ; ] ; 




Dry* +- 


louuuu^ri; 




al = 


[al; Rmt A 4 Rmt A 3 Rmt A 2 


Rmt A l; ] ; 


bl = 


[bl; eval (tmpmTx) ; ] ; 





These are the coeff. for the multi-torus eq. Rmt = variable 



% 
% 
% 

coeflS = pinv(al)*bl; 
disp( 'coef 15. . . 1 ) ; 



Rmt = 1200*rl; 
Rt = 20*rl; 

al = [Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 

bl = [eval (tmpmTx) - coef 15 (1) *Rmt A 4 - coef 15 (2) *Rmt A 3 - coef 15 (3) *Rmt A 2 coeflS (4 ) *Rmt / 
l;l; 

Rt = 40*rl; 

al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 

bl = [bl; eval (tmpmTx) - coef 15 (1) *Rmt A 4 - coef 15 (2) *Rmt A 3 - coef 15 (3) *Rmt A 2 coeflS (4) 
*Rmt A l;] ; 
Rt = 80*rl; 

al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 

bl = [bl; eval (tmpmTx) - coef 15 ( 1 ) *Rmt A 4 - coef 15 (2) *Rmt A 3 - coef 15 (3) *Rmt A 2 coeflS (4) 
*Rmt A l;] ; 
Rt = 120*rl; 

il = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 

ll = [bl,:.. eval (tmpmTx) - coerf!5 (1) *Rmt A 4 - coef 15 (2) *Rmt A 3 - coef 15 (3) *Rmt A 2 coeflS (4) 
'*Rmt A l;] ; 
%; 

%; These are the- coeff\ for the multi-torus eq. Rmt = 1200, Rt = variable 
%; 

coefl6 = pinv(al)*bl; 
disp ( 'coef 16. . . ' ) ; - 



%; Here is the final multi-torus equation (fox x here) 
%; 

syms multitoruslx multitorusly multitoruslz; 
coeflS = eval (coef 15) ; 
coef 16 = eval (coef 16) ; 

multitoruslx = coef 15 (1) *Rmt A 4 + coef 15 (2 ) *Rmt A 3 + coef 15 (3 ) *Rmt A 2 + coef 15 (4 ) *Rmt A l; 
multitoruslx = multitoruslx + coef 16 (1) *Rt A 4 + coef 16 (2 ) *Rt A 3 + coef 16 (3) *Rt A 2 + coefl6(4) 
* Rt A 1 ; 

disp ( 1 finished computation of multitoruslx...'); 
fprintf (f id, 1 \n\nmultitoruslx = Sum of: \n'); 
for il = 1:4 

fprintf (fid, 1 %20 . 12f 1 , coef 15 (il) ) ; 

fprintf (fid, ' Rmt A %li ' , il) ; 

f print f (fid, ' \n' ) ; 
end 

for il = 1:4 

fprintf (fid, ' %20 . 12f ' , coef 16 (il ) ) ; 
fprintf (fid, 1 Rt A %li',il); 
fprintf (fid, '\n' ) ; 
"end 

fprintf (fid, 1 \n\nTest of multitoruslx'); 
Rt = 20; 
Rmt = 100; 



f print f (fid, 1 \nRt = %20 . 12f 1 , Rt ) ; 
f print f (fid, 1 \nRmt = %20.12f f , Rmt) ; 

f print f (fid, ' \nmultitoruslx = %20 . 12f ' , multitoruslx) ; 
f print f (fid, '\n' ) ; 



Otxt = Rmt*cos (Thmt ) ; 
Otyt = Rmt*sin (Thmt ) ; 
Otzt = 0.0; 
Thtt = 0.0; 
nT = 16; 



tmpmTy = 0; 
for il = l:nT 

Thmt = il*2*3.14159/nT; 

tmpmTy = tmpmTy + toruseqly; 
end 



Rmt = 1000*rl; 

al = [Rmt A 4 Rmt A 3 Rmt A 2 Rmt A l; ] ; 
bl = [eval (tmpmTy) ; ] ; 
Rmt - 2000*rl; 

al = [al; Rmt A 4 Rmt A 3 Rmt A 2 Rmt A l;]; 
bl = [bl; eval (tmpmTy) ; ] ; 
Rmt = 10000*rl; 

al = [al; Rmt A 4 Rmt A 3 Rmt A 2 Rmt A l;]; 
bl = [bl; eval (tmpmTy) ;] ; 
Rmt = 160000*rl; 

al = [al; Rmt A 4 Rmt A 3 Rmt A 2 Rmt A l;]; 
bl = [bl; eval (tmpmTy) ;] ; 



% 



These are the coeff. for the multi-torus eq. Rmt = variable 



coef!5 = pinv(al)*bl; 
disp ( 1 coef 15 . . . f ) ; 

Rmt = 1200*rl; 
Rt = 20*rl; 

al = [Rt A 4 Rt A 3 Rt A 2 Rt A l;]; 

bl = [eval (tmpmTy) - coef 15 ( 1) *Rmt"4 - coef 15 (2 ) *Rmt A 3 - coef 15 (3) *Rmt A 2 coef 15 (4 ) *Rmt A 

l;]; 

Rt = 40*rl; 

al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 

bl = [bl; eval (tmpmTy) - coef 15 ( 1 ) *Rmt A 4 - coef 15 (2) *Rmt A 3 - coef 15 (3) *Rmt A 2 coef 15 (4) 
*Rmt A l; ] ; 
Rt = 80*rl; 

al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 

bl = [bl; eval (tmpmTy) - coef 15 (1) *Rmt A 4 - coef 15 (2) *Rmt A 3 - coef 15 (3) *Rmt A 2 coefl5(4) 
*Rmt A l; ] ; 
Rt = 120*rl; 

al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 

bl = [bl; eval (tmpmTy) - coef 15 (1) *Rmt A 4 - coef 15 (2 ) *Rmt A 3 - coef 15 (3) *Rmt A 2 coef 15 (4) 

*Rmt A l;] ; 

%; 

%; These are the coeff. for the multi-torus eq. Rmt = 1200, Rt = variable 
%; 

coefl6 = pinv(al)*bl; 
disp ( 1 coef 16. . . 1 ) ; 



%; Here is the final multi-torus equation 
%; 

syms Rt Rmt; 

multitorusly = coef 15 ( 1 ) *Rmt A 4 + coef 15 (2 ) *Rmt A 3 + coef 15 (3) *Rmt A 2 + coef 15 (4 ) *Rmt A l; 



multitor&sly = nrultitorusly + coef 16 (1) *Rt A 4 + coef 16 (2 ) *Rt A 3 + coef 16 (3) *Rt A 2 + coefl6(4) 
*Rt A l; 

disp( 'finished computation of multitorusly...'); 
fprintf (fid, 1 \n\nmultitorusly = Sum of: \n'); 
or il = 1:4 

fprintf (fid, ■ %20. 12f ? , coef 15 (il) ) ; 

fprintflfid, Rmt A %li f , il) ; 

fprintf (fid, '\n') ; 
end 

for il = 1:4 

fprintf (fid, f %20. 12f 1 , coef 16 (il) ) ; 

fprintf (fid, 1 Rt A %li',il); 

fprintf (fid, '\n f ) ; 
end 

fprintf (fid, 1 \n\nTest of multitorusly 1 ); 
Rt = 20; 
Rmt = 100; 

fprintf (fid, f \nRt = %20. 12f 1 , Rt) ; 
fprintf (fid, 1 \nRmt = %20.12f ' ,Rmt) ; 

fprintf (fid, 1 \nmultitorusly = %20. 12f f , eval (multitorusly) ) ; 
fprintf (fid, *\n' ) ; 



Otx = Rmt * cos (Thmt) ; 
Oty = Rmt * sin ( Thmt ) ; 
Otz = 0.0; 
nT = 16; 

tmpmTz = 0; 
for il = l:nT 

Thmt = il*2*3.14159/nT; 

tmpmTz = tmpmTz + toruseqlz; 
end 

Rmt = 1000*rl; 

al = [RmfM Rmt A 3 Rmt A 2 Hiatal; ]; 
bl = [eval (tmpmTz) ;] ; 
Rmt = 2000*rl; 

al = [al^^Mrt** Rmt A 3 Rmt A '2 Rmt A l;]; 
bl = [bl; eval (tmpmTz) ;] ; 
Rmt = 10000*rl; 

al = [alr ii HKtr*4 Rmt A 3 Rmt A 2 Rmt A l;]; 
bl = [bl; eval (tmpmTz) ;] ; 
Rmt = 160000*rl; 

al = [al; Rmt A 4 Rmt A 3 Rmt A 2 Rmt A l;]; 
bl = [bl; eval (tmpmTz) ;] ; 

%; 

%; These are the coeff. for the multi-torus eq. Rmt = variable 
%; 

coef 15 = pinv(al)*bl; 
disp( f coef 15. . . 1 ) ; 

Rmt = 1200*rl; 
Rt = 20*rl; 

al = [ Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 

bl = [eval (tmpmTz) - coef 15 ( 1 ) *Rmt A 4 - coef 15 (2) *Rmt A 3 - coef 15 (3) *Rmt A 2 coef 15 (4 ) *Rmt' 
l;]; 

Rt = 40*rl; 

al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 
1 = [bl; eval (tmpmTz) - coef 15 (1) *Rmt A 4 - coef 15 (2 ) *Rmt A 3 - coef 15 (3) *Rmt A 2 coef 15 (4) 
Rmt A l; ] ; 
t = 80*rl; 
al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 

bl = [bl; eval (tmpmTz) - coef 15 (1) *Rmt A 4 - coef 15 (2 ) *Rmt A 3 - coef 15 (3) *Rmt A 2 coef 15 (4) 
*Rmt A l;]; ^ 



Rt = 120*rl; 

al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A l;]; 

bl = [bl; eval(tmpmTz) - coef 15 (1) *Rmt A 4 - coef 15 (2) *Rmt A 3 - coef 15 (3) *Rmt A 2 coefl5(4) 



*Rmt A l;] ; 



These are the coeff. for the multi-torus eq. Rmt = 1200, Rt = variable 



coefl6 = pinv(al)*bl; 
disp( 'coef 16. . . 1 ) ; 



%; 

%; Here is the final multi-torus equation 
%; 

multitoruslz = coef 15 ( 1 ) *Rmt A 4 + coef 15 (2) *Rmt A 3 + coef 15 (3) *Rmt A 2 + coef 15 ( 4 ) *Rmt A l; 
multitoruslz = multitoruslz + coef 16 (1) *Rt A 4 + coef 16 {2 ) *Rt A 3 + coef 16 (3) *Rt A 2 + coef 16 (4) 
*Rt A l; 

disp( 'finished computation of multitoruslz...'); 
fprintf (fid, ' \n\nmultitoruslz = Sum of: \n ! ); 
for il = 1:4 

fprintf (fid, ' %20 . 12f ' , coef 15 (il ) ) ; 

fprintf (fid, ' Rmt A %li ' , il ) ; 

fprintf (fid, ■ \n f ) ; 
end 

for il = 1:4 

fprintf (fid, ' %20 . 12f ' , coef 16 (il) ) ; 

fprintftfid, '" Rt A %li ' , il ) r 

fprintf (fid, ' \n' ) ; 
end 

fprintf (fid, ' \n\nTest of multitoruslz'); 
Rt = 20; 
*mt = 100; 

fprintf (fid, ' \nRt = %20 . 12f ' , Rt ) ; 

fprintf (fid, ' \nRmt = %20 . 12f ' , Rmt ) ; 

fprintf (fid, ' \nmultitoruslz = %20 . 12f , multitoruslz) ; 
fprintf (fid, 1 \rr' ) ; 



compute the plasma-torus equation 

- plasma toroid equation is f (prl,pthl,pth2,pr2) 

- inner ring x = r*cos(pthl) y = r*sin(pthl) 

- outer branch is 1st toroid equation 

- vector inward is from outer branch point directed back to ring 

- final toroid equation is projection of Bx,By,Bz onto vector in 

- cos thbp = Bx*vix + circlely*viy + circlelz*viz/ | | B || | | vi | | 
%; - proj . = | | B | | cos thbp 

%; - ...the multi-torus crossed with the plasma-torus 

syms plrad plrad2 pthl plxl plyl plzl projl proj2 proj3 proj4 ptplxl ptplyl ptplzl; 

syms veclx vecly veclz; 

Otxt = Otx; 

Otyt = Oty; 

Otzt = Otz; 

Thtt = Tht; 

Thtnt = Thtn; 



pthl = 3.14159/2; 
plxl = plrad*cos (pthl) ; 
lyl = plrad*sin(pthl) ; 
lzl = 0.0; 

th2 - pthl/ (2*3. 14159*plrad) ; 
projl = plrad2*cos (pth2) ; 
proj2 = (plrad + proj 1) *cos (pthl) 
proj3 = plrad2*sin (pth2) ; 



proj4 = grtJji^sri-n (pthl) ; 

ptplxl = proj2; 

ptplyl = proj4; 

ptplzl = groj3;~ 
f^reclx = ptplxl - plxl; 
^K^ecly = ptplyl - plyl; 

veclz = ptplzl - plzl; 

plrad = Rmt; 

plrad2 = Rtt; 



syms disB disVec thbp projMt; 

disB = sqrt (multitoruslx A 2 + multitorusly A 2 + multitoruslz A 2) ; 
disVec = sqrt(veclx A 2 + vecly A 2 + veclz A 2); 

thbp = acos ( (multitoruslx*veclx + multitorusly*vecly + multitoruslz*veclz) / (disB*disVec) ) ; 
projMt = disB*cos (thbp) ; 

disp ( 1 finished computation of projMt... 1 ); 

fprintf (fid, 1 \n\ndisB = f ); 

fprintf (fid, 1 \ndisB - %20 . 12f » , eval (disB) ) ; 

fprintf (fid, f \ndisVec = f ) ; 

fprintf (fid, char (disVec) ) ; 

fprintf (fid, 1 \nthbp = '); 

fprintf (fid, char (thbp) ) ; 

fprintf (fid, 1 \nprojMt = '); 

fprintf (fid, char (projMt) ) ; 

fprintf (fid, ? \n f ) ; 



compute the minimum of the plasma-torus equation 

- compute the derivative of the plasma-torus equation 

- ...derivative is found with matlab 

- set the minimum to the necessary pressure and temperature for 

- ... fusion 

- cc ff lg Tule ' of f^re^-rc^rtrs of the derivative 

- ...roots are found with matlab 

- find the root that produces the minimum 

- . . Iffy goifig down" the * T rl^st 

f (Rt,Rmt, ptplxl, ptplyl', ptplzl) matrix (4 terms per variable) 
find t¥¥^^^ivat^e- 
find roots of first derivative 
syms coefl7 coefl8 coefl9 coef20 coef21; 
Otxt = Otx; 
Otyt = Oty; 
Otzt = Otz; 



Rt = 20; 
Rtt = Rt; 



al = [Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 

bl = [eval (projMt) ;] ; 

Rt = 40; 

al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 

bl = [bl; eval (projMt) ;] ; 

Rt = 80; 

k al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 

1 = [bl; eval (projMt) ;] ; 

.t = 160; 

al = [al; Rt A 4 Rt A 3 Rt A 2 Rt A 1 ; ] ; 

bl = [bl; eval (projMt) ;] ; 
coefl7 = pinv(al)*bl; 



disp( f coefl7. . . 1 ) ; 

Rt = 20; 
Rmt = 100; 

• 1 = [Rmt A 4 Rmt A 3 Rmt A 2 Rmt A l;]; 
1 = [eval(projMt) - coef 17 (1) *Rt A 4 - coef 17 (2) *Rt A 3 - coef 17 (3) *Rt A 2 - coef 17 (4 ) *Rt A l; ] ; 
Rmt = 200; 

al = [al; Rmt A 4 Rmt A 3 Rmt A 2 Rmt A l;]; 

bl = [bl; eval(projMt) - coef 17 (1) *Rt A 4 - coef 17 (2) *Rt A 3 - coef 17 (3) *Rt A 2 - coef 17 (4 ) *Rt A 
l;]; 

Rmt = 400; 

al = [al; Rmt A 4 Rmt A 3 Rmt A 2 Rmt A l;]; 

bl = [bl,v eval(projMt) - coef 17 (1) *Rt A 4 - coef 17 (2 ) *Rt A 3 - coef 17 (3) *Rt A 2 - coef 17 ( 4 ) *Rt A 
l;]; 

Rmt = 800; 

al = [al; Rmt A 4 Rmt A 3 Rmt A 2 Rmt A l;]; 

bl = [bl; eval(projMt) - coef 17 (1) *Rt A 4 - coef 17 (2) *Rt A 3 - coef 17 (3) *Rt A 2 - coef 17 (4) *Rt A 
l;]; 

coef 18 = pinv(al)*bl; 
disp( 'coef 18. . . 1 ) ; 



syms plasmatorusl plasmatorusdif fl; 

syms plasmatorusdif f2Rt plasmatorusdif f2Rmt; 

syms plasmatorusdif f3Rt plasmatorusdif f3Rmt; 

syms sizel size2 iil ii2 curmin testmin minRt minRmt; 

plasmatorusl = coef 17 (1) *Rt A 4 + coef 17 (2 ) *Rt A 3 + coef 17 (3) *Rt A 2 + coef 17 (4) *Rt A l; 
plasmatorusl = plasmatorusl + coef 18 (1) *Rmt A 4 + coef 18 (2) *Rmt A 3 + coef 18 (3) *Rmt A 2 + coefl8 
(4)*Rmt A l; 

disp ( 1 finished computation of plasmatorusl . . . 1 ) ; 
fprintf (fid, 1 \n\nplasmatorusl = Sum of: \n ! ); 
^or il = 1:4 

fprintf (fid, ' %20 . 12 f f , coef 17 (il) ) ; 

fprintf (fid, 1 Rt A %li',il); 

fprintf (fid, 1 \n' ) ; 
end 

for il = 1:4 

fprintf (fid, f %20 . 12f 1 , coef 18 (il) ) ; 

fprintf (fid, 1 Rmt A %li 1 , il ) ; 

fprintf (fid, f \n f ) ; 
end 



plasmatorusdiff2Rt = [coefl7(l)*4 coefl7(2)*3 coef 17 (3) *2 coef 17 (4) *1] ; 
plasmatorusdif f2Rmt = [coefl8(l)*4 coefl8(2)*3 coefl8(3)*2 coef 18 (4) *1] ; 

plasmatorusdif f3Rt = roots (plasmatorusdif f2Rt) ; 
plasmatorusdif f3Rmt = roots (plasmatorusdif f2Rmt) ; 
fprintf (fid, 1 \n\nplasmatorusdiff 3 = Sum of: \n f ); 
for il == 1:3 

fprintf (fid, ! %20 . 12f 1 , plasmatorusdif f3Rt (il) ) ; 

fprintf (fid, 1 plasmatorusdif f3Rt A %li 1 , il ) ; 

fprintf (fid, 1 \n f ) ; 
end 

for il = 1:3 

fprintf (fid, 1 %20 . 12f 1 , plasmatorusdif f3Rmt (il ) ) ; 
. fprintf (fid, 1 plasmatorusdif f 3Rmt A %li 1 , il) ; 
I fprintf (fid, 'Xn 1 ); 
^nd 



sizel = size (plasmatorusdif f3Rt, 1) ; 



size2 = size (plasmatorusdif f 3Rmt, 1) ; 



syms plamatl Rtx Rtmx ptplxlx ptplylx ptplzlx; 

plasmatl = coef 17 ( 1) *Rtx A 4 + coef 17 (2) *Rtx A 3 + coef 17 (3) *Rtx A 2 + coef 17 (4 ) *Rtx A l; 
lasmatl = plasmatl + coef 18 (1) *Rtmx A 4 + coef 18 (2) *Rtmx A 3 + coef 18 (3) *Rtmx A 2 + coefl8(4) 
Rtmx A l; 



curmin = 9999999; 
for iil = Irsizel 
for ii2 = l:size2 

Rtx = plasmatorusdif f3Rt (iil) ; 

Rtmx = plasmatorusdif f3Rmt (ii2) ; 

testmin = eval (plasmatl) ; 

if or ((testmin < curmin), (curmin == 9999999)) 
curmin = testmin; 
minRt = Rtx; 
minRtm = Rtmx; 
end 
end 
end 

disp ( 1 finished computation of minimum of plasmatorus eq. . . 1 ) ; 
fprintf (f id, 1 \n\nMinimum of Plasmatorus Eq Occurs At:'); 
fprintf (fid, 1 \nminRt = %20. 12f 1 , minRt) ; 
fprintf (fid, ' \nminRtm = %20. 12f f , minRtm) ; 
fprintf (fid, f \n») ; 

syms Rtq Rtmq plasmatorusm; 
Rt = minRt; 
Rmt = minRtm; 

plasmatorusm = coef 17 ( 1) *Rtq A 4 + coef 17 (2 ) *Rtq A 3 + coef 17 (3) *Rtq A 2 + coef 17 (4 ) *Rtq A l ; 
plasmatorusm = plasmatorusm + coef 18 (1) *Rtmq A 4 + coef 18 (2 ) *Rtmq A 3 + coef 18 (3) *Rtmq A 2 + 
goefl8(4)*Rtmq A l; 

Plisp ( 1 finished computation of polynomial eq. of min of plasmatorus eq. . . 1 ) ; 
fprintf (fid, f \n\nplasmatorusm = Sum of: \n f ); 
for il = 1:4 

fprintf (fid, 1 %20 . 12f ' , coef 17 (il) ) ; 

fprintf (fid, 1 Rt A %li',il); 

fprintf (fid, 1 \n ! ) ; 
end 

for il = 1:4 

fprintf (fid, 1 %20 . 12f 1 , coef 18 (il) ) ; 

fprintf (fid, 1 Rmt A %li 1 , il) ; 

fprintf (fid, f \n f ) ; 
end 



%; 

%; general equation of surface of a torus with radiuses Rt and Rtm 
%; 

syms Surfacel Rtmm Rtmmm; 

Surfacel = 4*3 . 14159*3 . 14159*Rtmm*Rtmmm; 

syms plasmatm Rtqq Rtmqq ptplxlqq ptplylqq ptplzlqq; 

plasmatm = coef 17 ( 1) *Rtqq A 4 + coef 17 (2) *Rtqq A 3 + coef 17 ( 3 ) *Rtqq A 2 + coef 17 ( 4 ) *Rtqq A l ; 
plasmatm = plasmatm + coef 18 ( 1 ) *Rtmqq A 4 + coef 18 (2 ) *Rtmqq A 3 + coef 18 (3) *Rtmqq A 2 + coefl8 
(4)*Rtmqq A l; 

m, 

%; compute the current equation for ring sizes 

%; - polynomial I (10) Rl (10) R2 (10) R3 (10) Nl (10) N2 (10) 



compute a table of current, ring sizes, number of rings, etc. 

- for current (1 Amp to 100 Amp) (wait on this) 

- ring size 1 {1 to 5 meters) 

- ring size 2 (10 to 20 meters) 

- ring size 3 (40 to 80 meters) 

- # of rings 1 (10 to 20) 

- # of rings 2 (10 to 20) 

- # of rings 3 (10 to 20) 

some constants: 

mass of electron = 9.110 * 10 A -31 kg 
%; charge of electron = 1.602192 * 10 A -19 C 
%; solenoid constant (uO) - 4*pi*10 A -7 T A A -1 
%; solenoid equation: B = uO * n * I 
%; Gas constant = 8.314 J mol A -l K A -1 
%; Avogadro's # = 6.022 * 10 A 23 molecules /mole 
% 

syms melee qelec soluO gask avogad speedc vor; 

melee = 9. 110*10 A (-31 ) ; 

qelec = 1 . 602192*10 A (-19) ; 

soluO = 4*3.14159*10 A (-7) ; 

gask = 8.314; 

avogad = 6.022*10 A (23) ; 

speedc = 2 . 9979*10 A (8) ; 

vor = 100/gask; 

syms rsizel rsize2 rsize3 nringsl nrings2 minB fperA curl cur2; 

syms Efmt Ffmt vfmt nfmt Bpar Fpar Vpar Energyl Pressurel Temperaturel; 

%; compute and print: 

%; Efmt = speedc*plasmatm; 
%; Ffmt = qelec*Efmt; 
%; vfmt = (Ffmt/melec) *60; 
%; nfmt - 2*3.14159/0.05; 
Bpar - nfmt*cur2*solu0; 
Fpar = qelec*vfmt*Bpar; 
%; vpar = (Fpar/l)*l; 
%; Energyl = ( 1/2 ) *l*vpar A 2/vor; 
Temperaturel = gask*Energyl ; 

Pressurel = vor*gask*Temperaturel/ (4*3 . 14159*3 . 14159*rsize3*rsize3) ; 
f print f (fid, ' \n\nTable of Minimum B Field Due to Various Radius Sizes and Number of 
Toruses 1 ) ; 
f print f (fid, 1 \n ? ) ; 
for rsizel = 1:2 
for rsize2 = 10:10 
for rsize3 = 40:40 
for nringsl = 10:11 
for nrings2 = 10:11 
for nrings3 = 10:10 
for curl = 1:1 
for cur2 =1:1 
rl = rsizel; 
Rtqq = rsize2; 
Rtmqq = rsize3; 
nr = nringsl; 
nT = nrings2; 
Rtmm » Rtmqq; 
Rtmm = rsizel; 
Rtmmm = rsize2; 
minB = eval (plasmatm) ; 
fperA = minB*eval (Surfacel) ; 
Efmt = speedc*plasmatm; 
Ffmt = qelec*Efmt; 
vfmt = (Ffmt/melec) *60; 
nfmt = 2*3.14159/0.05; 
Bpar = nfmt*cur2*solu0; 
Fpar = qelec*vfmt*Bpar; 
vpar = (Fpar/1) *1; 



Energyl = (1/2) *l*vpar A 2/vor; 
Temperaturel = gask*Energyl; 

Pressurel = vor*gask*Temperaturel/ (4*3. 14159*3. 14159*rsize3*rsize3) ; 

m, 

%; required fperA is a constant from textbook 

%; 

disp ( ' inner loop ...'); 

rsizel 

rsize2 

rsize3 

nringsl 

nrings2 

fprintf (fid, ' \n' ) ; 

f print f (fid, 'Ring Size 1 = '); 

fprintf (fid, 1 %20 . 12f ' , rsizel ) ; 

fprintf (fid, '\n' ) ; 

fprintf (fid, 'Ring Size 2 = ! ); 

fprintf (fid, ' %20 . 12f ' , rsize2 ) ; 

fprintf (fid, '\n f ) ; 

fprintf (fid, 'Ring Size 3 = '); 

fprintf (fid, 1 %20 . 12f 1 , rsize3) ; 

fprintf (fid, '\n' ) ; 

fprintf (fid, '# of Rings 1 = ! ); 

fprintf (fid, 1 %4i 1 , nringsl) ; 

fprintf (fid, f \n f ) ; 

fprintf (fid, f # of Rings 2 = f ); 

fprintf (fid, 1 %4i 1 , nrings2 ) ; 

fprintf (fid, ? \n' ) ; 

fprintf (fid, 'minB = '); 

fprintf (fid, 1 %20 . 12f 1 , minB) ; 

fprintf (fid, 'Xn 1 ) ; 

fprintf (fid, 1 fperA = '); 

fprintf (fid, f %20. 12f ' , fperA) ; 

fprintf (fid, f \n') ; 

fprintf (fid, 1 \n f ) ; 

fprintf (fid, 'curl = '); 

fprintf (fid, ' %20 . 12f ' , curl ) ; 

fprintf (fid, ' \n' ) ; 

fprintf (fid, 'cur2 « '); 

fprintf (fid, '%20.12f ' , cur2) ; 

fprintf (fid, ? \n' ) ; 

fprintf (fid, 'Temperaturel = '); 

fprintf (fid, ' %20 . 12f ', eval (Temperaturel ) ) ; 

fprintf (fid, ' \n' ) ; 

fprintf (fid, ' Pressurel = '); 

fprintf (fid, ' %20 . 12f ', eval (Pressurel) ) ; 

fprintf (fid, ' \n' ) ; 



print both of them for comparison 



end 
end 
end 
end 
end 
end 
end 
^nd 

disp (' finished program, 
st = fclose (f id) ; 



